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Set Books 


G. D. Smith, Numerical Solution of Partial Differential Equations (Oxford. 1971). 
Н. F. Weinberger, A First Course in Partial Differential Equations (Xerox, 1965). 


It is essential to have these books; the course is based on them and will not make 
sense without them. They are referred to in the text as Sand W respectively. 


Unit 15 includes a short reading passage from 5, but is not otherwise based on either 
set book. 


Conventions 


Before working through this text make sure you have read А Guide to the Course: 
Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 


Unit M100 13, Integration II for the Mathematics Foundation Course, 
Unit M201 23, The Wave Equation for the Linear Mathematics Course. 
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15.0 INTRODUCTION 


In the three previous numerical units we have looked at some simple partial differential 
equations and developed finite-difference schemes to solve them. We have been 
interested primarily in the conditions under which our numerical methods can be 
used with confidence. This has led us to discuss stability and convergence in the 
solution .of initial value problems and convergence of iterative methods in various 
contexts. 


Inthis unit we shall look at some of the problems we have neglected, dealing exclusively 
with parabolic equations. 


We begin the unit by attempting to improve the accuracy.of finite-difference schemes 
applied to parabolic equations in one space dimension. This investigation leads us 
naturally into a discussion of three-level schemes. It so happens that three-level 
schemes are useful in overcoming some of the difficulties associated with the numerical 
solution of nonlinear partial differential equations and so we devote much of the 
discussion to' the treatment of a simple class of such equations. This section will 
illustrate the usefulness of finite-difference schemes in an area where it is extremely 
difficult to obtain analytical solutions. We shall also devote a short section to equations 
in two dimensions. 


The final section is a discussion of what happens when the data involve discontinuities. 
In previous work, we have been concerned mainly with problems in which the 
differential equations have no discontinuities in their coefficients. and in which the 
subsidiary conditions are differentiable. There are, of course, problems in which 
these conditions are not satisfied and we shall discuss briefly the modifications which 
must be made to our finite-difference schemes to deal with these situations. As one 
might expect, there is no universal panacea for such problems and we have to rely 
heavily on our knowledge of the analytical properties of partial differential equations 
to give us some insight into possible lines of approach in individual cases. Indeed, 
it is always important to consider analytical and numerical techniques in conjunction 
with each other, for they can interact to suggest suitable techniques of solution. 


151 PARABOLIC EQUATIONS IN ONE SPACE DIMENSION 


15.1.1 Three-Level Schemes 


In our earlier numerical work on parabolic equations, we concentrated on two-level 
schemes; that is, schemes which involve combinations of pivotal values on two time 
levels. We first introduced explicit schemes, which are easy to compute but give rise 
to problems of instability. Stable explicit methods require small step sizes in the 
t-direction and hence involve much computational effort. Implicit methods, on the 
other hand, give more involved computations at each step, but they allow us to use 
larger values of the step size because of good stability properties. One of the most 
important implicit methods is the Crank-Nicolson method which, as well as being 
unconditionally stable, has a small local truncation error. It should be clear that to 
be able to reduce the size of local truncation error of a scheme whilst retaining 
stability is of great value in the numerical solution of partial differential equations. 
This is because we can either obtain better accuracy in a given number of steps or 
get similar results with less effort by taking larger steps. 


In the case of the simple diffusion equation, 
aU _ eu 
a D 


> 


we saw in Unit 5, Initial Value Problems that replacement of the time derivative using 
the forward-difference formula 


where k is the mesh spacing in the r-direction, leads to the simple explicit scheme 


Шука = ш + шарр 20; + ti-a.) (1) 
where r = k/h?. If we were to replace the time derivative using the more accurate 
central-difference formula 

QU, ; VUL je 7 Hj 

t 2k 


then we would obtain the following explicit scheme: 


Шіл 2ғ(ш 1; - 2; ta) Шшу-і: (2) 


It is immediately clear that this is a three-level scheme since it involves terms along the 
Q — Ith, jth and (j + 1)th time levels. Our quest for higher accuracy has led us 
naturally to a three-level scheme and we must now see whether we have in fact pro- 


duced a Scheme with a smaller local truncation error. This is the subject of the 
following SAQ. 


SAQI 


Show that the simple three-level explicit scheme given by Equation (2) is locally 
more accurate than the two-level explicit scheme gi i ) 

; given by Equation (1) whi Г 
are applied to the diffusion equation. 2 шс; 


(Solution on p. 24.) 


Although we have produced a more accurate explicit scheme, it is of little use unless 
it turns out to be stable for some values of the mesh ratio r, and it is this aspect of 
three-level schemes which we shall now discuss. Rather than look at the explicit 


scheme on its own we shall include an а: i ic ri = 
nalysis which works al: i ici 
vi | so for implicit three: 


By analogy with the work in Section 8.4.1 of Unit 8, Stability you can see that a general 
three-level scheme can be written in matrix form as 


Аш, = Bu; + Cuj., (3) 
or, provided det А # 0, as 
ш; = A^! Buj + A7) Cuj. ,. 


(For simplicity in the analysis we are considering the case of homogeneous boundary 
conditions, that is. U(0, t) = U(1.r) = 0.) By putting v;,, = uj we can express this 
three-level scheme as a pair of two-level schemes, which can be written in partitioned 
matrix form as 


ui] [478 АТС 
bos Е 


where 1 is the unit matrix of the same order as the matrices А, B and C. We call the 
matrix P, given by 


АЗВ ACC 
P- И 
1 0 
the amplification matrix, and we know from Section 8.4.1 of Unit 8, which dealt 
with the matrix method for examining stability, that the scheme (4) is stable provided 


the spectral radius of the amplification matrix P is not greater than unity. Therefore, 
we now need to find the eigenvalues of P. If 4 is an eigenvalue of P then we may write 


A^ B  A^C][u .[u 

а 
where (и, v) is ап eigenvector. Expanding the equation, we obtain 

A7! Bu + A7! Cv = 2u, 

u = дү, 

Eliminating u between the equations, we are left with 

GI — АТВ — A7! Ov = 0, 
and the eigenvalues of P satisfy 

det?] — 44^! B — A7!C) = 0, 
whence 

det(24 — åB — C) = 0. 
It is this equation which we find the most useful for determining the eigenvalues of 
the matrix P. 
Example 


Determine the eigenvalues of the amplification matrix of the three-level explicit 
scheme given by Equation (2) for the diffusion equation with homogeneous boundary 
conditions. Hence investigate the stability of this scheme. 


Solution 
Writing the scheme in the form of Equation (3) we see that А = C = I and Bis given by 
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The eigenvalues of the amplification matrix are given by 
det(/?7 - 4B — 1) 20 
which, on dividing by 2 (assuming 2 # 0), gives 


1 


det [s - i =0. 
5 


Therefore. if j is an eigenvalue of B.then 2 satisfies 


The eigenvalues of B are given by 


«= ват (Z) з-1,2,....М-1 


where B has M — І rows and columns (see Appendix to Unit 8), and so 


22 + Agrsin? (| -1=0. 


2M, 
Hence, 
А sol sm 3 T 
АЁ = —4r sin? (| + J1 + 162 sin* (=| &m12,.M = 1, 
and since 4; < —1 for all s the scheme is unstable for all values of r. 


This example has shown us that the three-level explicit formula, whilst of the same 
local accuracy as the Crank- Nicolson scheme, and easier to use because it is explicit, 
fails because of instability. There are, however, successful formulas which are accurate, 
explicit and stable as we ask you to show in SAQ 2. 


SAQ 2 


Use the matrix method to investigate the stability of the explicit Du Fort and Frankel 
scheme, 


Шыжа T Uij-a Шар Муж — 8-і t Minny 


2k h? 
applied to the problem 
UPU Gone 0 
x a «x«l, г> 0, 


О(0, N = U(1.) = 0 г> 0. 
(Solution on р. 24.) 


In SAQ 7 of Unit 8 we showed that the local truncation error of the Du Fort and 
Frankel scheme applied to the simple diffusion equation is 0(42) + O(h?) + O(K?/h?), 
and that the scheme is consistent with the diffusion equation provided 


im - = 0. 

hk~0 h 

Since our last SAQ has shown the method is stable, it is also convergent (by Lax's 
Theorem), and we can now ask whether it is better than the Crank-Nicolson scheme. 
Since the local order of accuracy of the Crank- Nicolson method is O(h?) + O(k?) 
we cannot compare the two schemes directly because of the term O(k?/h?) in the 


local truncation error of the Du Fort and Frankel scheme. The following SAQs 
provide some insight into the manner of resolving this problem. 
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SAQ 3 
(a) Find the maximum order of accuracy of the Du Fort and Frankel scheme. 


HINT: Set k = O(h”) as h =з 0, where p is some real number, and plot the order 
of accuracy as a function of p. 


(b) Show that p — 2 is the most beneficial choice in (a). What does this result imply 
when choosing suitable values for the mesh ratio г = k/h?? 


(Solution on p. 26.) 


SAQ4 


By setting k = O(h”) as h -«0 in the formula for the local order of accuracy of the 
Crank-Nicolson scheme, determine the most suitable value for p. 


(Solution on p. 26.) 


The two previous SAQs have shown us that optimum use is made of the Du Fort 
and Frankel scheme if we choose К = O(h?) whereas the Crank-Nicolson scheme 
is best when used with k — O(h). For such optimum use both schemes have a local 
order of accuracy equal to O(h?). Therefore, if we consider M — 1 internal mesh 
points along each time level, so that h = 1/M, the step length k for the Du Fort and 
Frankel scheme will have to be of the order of 1/M?. Therefore, to step forward one 
unit in time with the Du Fort and Frankel scheme requires something like M? steps. 
A similar calculation with the Crank-Nicolson scheme shows that only about M 
steps are required for the same advance in time. 


In order to compare these two schemes fairly we need to find the amount of work 
required by each scheme to evaluate a solution, along some fixed time level, with 
comparable accuracy. We have just seen that more steps are required with the Du 
Fort and Frankel scheme than with the Crank-Nicolson scheme to reach a fixed 
time level with the same local order of accuracy. However, the Crank-Nicolson 
scheme is implicit and requires more computation at each step. To begin we need to 
estimate the amount of computation required by each scheme in advancing one 
step in time. Writing the Du Fort and Frankel scheme in the form 


a 
Шыжа = Mosa c LS (uei — js, Z Mig t Mima gh 
1+а 


where a = 2k/h?, we see that to evaluate Uu; j+ We need to perform one multiplication 
and four additions or subtractions. Therefore. the calculation of a complete level of 
values, assuming M — | internal mesh points along each time level, involves M — 1 
multiplications and 4M — 4 additions or subtractions. The Crank- Nicolson scheme 
given by 


TM ja H (2 + 27) ща ja, — uua jen = Tua, + (2 — Quay t ruis 


gives rise to a system of linear equations which can be written in the following form, 


assuming M — 1 internal mesh points along each time level: 
b, =c Uy jet d, 
та b; =n Ma jet dy 
—ам-› Әм-і —Cu-2 Uy 3j du- 
Чм -ај+а йм-1 
-а4м-і Әм-а d 
where b; = 2 + 2r. a; = c; = rand 


d; = rui, + (2 — 2ru j rui. 


These equations can be solved by the recurrence relation form of Gauss elimination. 
Referring to S: pages 21 and 22, we сап see that this method of solution will involve 
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5M — 9 multiplications or divisions and ЗМ — 6 additions or subtractions (notice 
that we need only evaluate the a;/z;.. , once). The calculation of the d; requires 2M — 2 


additions and 2M multiplications. (Note that. for example. 
d; = ruj- + (2 — 20m; + тщр 

and 
diaz = Tisi T (2 — 202; + ries 


the term ғи; ,.; occurs in both and need only be calculated once.) The amount of 
work required by the Crank—Nicolson scheme in stepping forward one level in time 


is therefore 7 — 9 multiplications or divisions and 5M — Sadditions or subtractions. 


In performing M? steps with the Du Fort and Frankel scheme (to step forward one 
unit in time), we shall require (M — 1)M? multiplications and 4(M — 1)M? additions. 
To reach the same time level with the Crank-Nicolson formula it takes of the order 
of M steps with a total of (7M — 9)M multiplications or divisions and (5M — 8)M 
additions or subtractions. 


To compare these figures we can assume that a modern computer takes the same 
time to compute a multiplication as it does a division and that additions take the 
same time as subtractions. Furthermore. for our present purposes. we shall assume 
that it takes ten times as long to compute a multiplication as it does an addition. 
With these assumptions we can easily deduce that for M > 5 (approximately) the 
Crank-Nicolson scheme is more efficient than the Du Fort and Frankel scheme. 


SAQ 5 
Compare the Crank—Nicolson and Du Fort and Frankel schemes for the solution of 
ðU QU 
r^ ЭХЕ 0«x«l, t>0, 
U(x, 0) = sin zx O<x<l, 
U(0.) -U(lL0)-20 і>0 


by obtaining values of u along the level t = 1. Use several values of h with the cor- 
responding optimum values of k obtained in SAQs 3 and 4. 


Show that the choice of k = h in the Du Fort and Frankel scheme results in a large 
error. How do you account for this? 


(Solution on p. 27.) 


It is also possible to construct implicit schemes. which do not have the consistency 
problem of the Du Fort and Frankel scheme (SAQ 7 of Unit 8). For example, we can 
construct a three-level implicit scheme in the manner in which the two-level Crank- 
Nicolson method was created. That is, we replace the second-order derivative 
0^ U/éx* by the mean of its central difference representations on the three levels 
j+ 1,jandj — 1. The three-level, Crank-Nicolson type. implicit replacement of 


eU у 
et ôx? 
is 
Hiper Маја _ 1 OU: зъл + біш; + бішу-а 
2k 3 h? 
or 
а- Збіш jay = $róu,; ж(- %б2ш;-1. (5) 


This scheme is stable for all values of r and has a local order of 


On) + O(K2). accuracy equal to 


When a three-level scheme is used, initial data must be provided on levels j = 0 and 
із 1 іл order to start the computation. Since a properly posed initial-boundary 


10 
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value problem of parabolic type provides initial data along г = 0 only, we need to 
generate values along t = k (i.e. j = 1) by some other means before we can use the 
three-level scheme. We could use a two-level difference formula of comparable 
accuracy to the three-level scheme or alternatively look for some analytical technique. 


SAQ 6 
Which scheme would you prefer to use to solve the simple diffusion equation 
20 CU. 
ôr аҳ? ` 


(a) the two-level Crank-Nicolson scheme. 
or 
(b) the three-level implicit scheme given by Equation (5)? 


(Solution on p. 28.) 


SAQ7 


A three-level implicit scheme for the diffusion equation, which has the same local 
order of accuracy as the Crank-Nicolson scheme and is unconditionally stable, is 
given by 


З [шуу — uu; l[ui; щу ШЕ 
1 =: Br ud | ip бына: 


Draw its molecular diagram. Show that this scheme involves less computational 
effort than the Crank-Nicolson scheme. 


(Solution on p. 28.) 


It is possible to construct (implicit) three-level schemes of higher order of accuracy 
which are always stable. However. such schemes involve more mesh points (typically 
up to nine points) and hence require more computational effort. 


At first sight it would appear that threc-level schemes have little more to offer us in 
the way of efficiency than two-level schemes. particularly since we have (usually) to 
use a two-level scheme to start them off. This is certainly true with linear equations. 
However, as we shall see in Section 15.1.4, three-level schemes can be of use in the 
solution of nonlinear equations. Before we go on to this topic we shall point out 
some of the difficulties associated with the numerical solution of nonlinear partial 
differential equations. 


15.1.2 Nonlinear Equations 


We have previously motivated our study of finite-difference methods by saying that 
at the present time it is usually only by numerical methods that we can effectively 
tackle nonlinear partial differential equations. The aim of this section is to present 
an introduction to this topic. We shall see that three-level schemes which we discussed 
in the previous section are of importance here. However, unlike the linear case, there 
is no unified theory for the solution of nonlinear partial differential equations. It 
turns out, however, that many techniques developed for linear equations apply 
equally well to nonlinear cases. 


Here we shall deal only with particular initial-boundary value problems of “parabolic” 
type given by 

ðU 20 ги 

rre убу 0<х<1, 0<:<Т (1) 
where f is some suitably well behaved function, together with (possibly) nonlinear 
boundary conditions.t When this problem is subject to smooth initial and boundary 
conditions it is properly posed in the region 0 < x < 1, 0 < t < T provided that 


of 
502 >а (2) 


Гог some constant а > 0. (The notation df /ĝU ,, represents the first partial derivative 
of the function f with respect to its fifth argument, 9? U/0x?.) 


Notice that when the right-hand side is just 9? U/8x? so that Equation (1) reduces to 
the simple diffusion equation, we have 


af _ 
QU. 


and Equation (2) is clearly satisfied. 


1 


Many of the numerical methods and techniques of proof for linear equations with 
constant coefficients carry over to nonlinear equations. However, in general, both the 
numerical processes and questions of stability and convergence are more complicated. 
The latter questions are largely unsolved and we generally replace "stability" by 
"local stability" on the assumption that the true solution of the differential equation 
varies little over small regions. We shall illustrate this point by investigating the 
stability in a simple example. 


For demonstration purposes we consider the nonlinear equation 


20 _ әу") 
зг әх? 9) 


where m is some integer greater than 1. The simplest explicit scheme for this equation 
is given by 


m 
Шужа Шыр Mea — 2ш + LESWA 


к h (4) 


The investigation of stability involves finding out how any errors in the pivotal values 
along one time level propagate to the next level. Suppose now that we perform the 
computation of Equation (4) with the values už; instead of u 
duced an error E; ; into u; ; given by ; 


ij» Where we have intro- 


= ut, — 
E; — ej tij 


t The analytical treatment of nonlinear partial differential equations does not form part of this course; 
we shall not pursue further the properties that the function f should possess to make the problem tractable. 


We feel obliged to point out, however, that considerable difficulties can arise in dealing with nonlinear 
differential equations and great care must be exercised. 
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That is, we evaluate an approximation to t; j+1 by the formula 


Uriel = ub _ (иё 1)" = 20и)" + (ut. Йй". 


k h? (5) 


Subtracting Equation (4) from Equation (5) and putting r = k/h?, we obtain 
(ші = шул) — (ut; — ш) = {иб д” = iay] 200)" — ш") 
+ uu}. (9 


Now, if we assume that E, ; is sufficiently small (i.e. E?; =~ 0 for m > 2), we can drop 
terms of second order and above in the binomial expansion. and obtain 


(ау — и, (щу + Е)" – ит ит + mE, ущ" up = mE, jury. 


We сап obtain similar expressions for (uf, , j)" — u; j and (uf, J" — ujy 80 
that Equation (6) becomes 
Eiji — Еултбшу Epa y 20 Ej + ШЕ.) (7) 


j 


ъ А 
The assumption that и varies little over a small region means that u;, 1j% uj% 4-1. 
and so Equation (7) gives 


Eijs = EQ теш (Ei, - 2Е;; + EL, 


We now apply von Neumann's method for investigating stability to this equation; 
we substitute into the equation the identity 

Eng = eil phia. (i NES 1) 
which yields the conclusion, eventually, that the process is stable if 

mr! <4, (8) 
provided that ий! > 0. (The scheme is unconditionally unstable if иг! < 0) This 
elementary analysis suggests that, for nonlinear problems, stability depends not only 
on the mesh ratio r but also upon the solution being obtained. (Recall that, in the case 
of the (linear) one-dimensional heat equation the condition r < } is necessary for 
stability, as we showed in Section 8.4 of Unit 8.) 


In practice, when solving nonlinear partial differential equations, we keep a constant 
check on stability (if the difference scheme is not unconditionally stable) by making 
tests, using a condition such as inequality (8), and either stop when the condition is 
not satisfied or automatically alter k to restore stability. 


Returning to Equation (1), a satisfactory explicit finite-difference replacement can 
be written as. . 


{зт ц; 75-і; 
шуу = uy ЖЕЛІК jk ц, з =y 


1 
© біш: 
This approximation is very simple to use, but suffers from the disadvantages that the 
mesh ratio r is strictly limited for stability. 


As in the linear case, the stability restriction can usually be removed by turning to 


implicit schemes of the Crank- Nicolson type. For Equation (1), the Crank-Nicolson 
equivalent scheme is 


Uje ij 1 
_ т = PUE + fj 9) 


with obvious notation. 


Unfortunately, depending upon the form of f, the algebraic problem of determining 
щӊу+ from Equation (9) may become quite complicated, since, if f is nonlinear, 
the algebraic equations at each time level are usually nonlinear. We may still be able 
to solve these equations however, as we shall illustrate in the following section. 
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5408 


Making the same assumptions as in the text show that the Crank-Nicolson type 
formula (9) is unconditionally stable locally when applied to Equation (3). (Use the 
von Neumann method and assume up! z0) 


(Solution on p. 29.) 


15.1.3 Iterative Solution of Nonlinear Algebraic Equations 


In most cases it is not possible to give an expression in closed form for the solutions 
of systems of nonlinear algebraic equations. Even when such expressions exist they 
may not be very useful for numerical purposes. As is frequently the case in numerical 
analysis, iterative methods which lead to approximations for the solutions are the 
best tools for this problem. There exist a large number of different methods, some 
better than others, but no one method is best in all cases. We shall briefly look at 
one such method which is often used, namely Newton's method. 


In Unit M100 14, Sequences and Limits П Newton’s method for solving the single 
equation 


f(x) =0 

was presented. It is the iterative scheme given by 
(nouo SO") = 
xt =a” т) n = 0,1,2,.... 


If the initial iterate x! is sufficiently close to a root о of f (x) = 0 then the successive 
iterates x"), x), |, converge to 0. 


Let us now look at the case of two equations in two unknowns given by 
Лб, ха) = 0 
Аха, ха) = 0. 


вов that we have initial approximations x) and xf? to the roots o, and ж, such 
that 


a, = x + dx, 


аз = xP + бх,. 


Then first-order Taylor approximations (Unit M20/ 14, Bilinear Forms) give the 
relations 


Ла, 0) = fio? + бху, x) + ôx) 
д 
= f xl), x!) + 2h um. x) 5x, + Л до, x9)óx, 
1 Ха 
апа 
Says a) = fa(xtD + óx,, x + 5x2) 
9f. a: 
~ f(x, 50) 2 (40) V fo 
= fX, xP) + ax, x, + ax xff)àx, . 


ees these expressions equal to zero, for a solution, we can write them in matrix 
orm as 


әл of 

0 0: қ 

femp е СЕС 
+ 


(01, xt Da Da = 
Лбх", x2”) ox; ax, 5x2 0 


co 
>= 
c 
> 
I 


where the derivatives are evaluated at (x. x), We see th 


(óx,. 5x5) is given by 


бх\ 


бх; 


where 


ox, 


wg 
БЕЛЕСІ 


an 
6x; 
f, 


Ox 


Лу, xP) 


ЛО, x) 
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at the correction vector 


The notation /7!\„ х,у means that the matrix J^! is evaluated using the values x, 
and ху. The iterative scheme is easily seen to be given by 


xr] [an лоб, ade) 


mE 
- = JT eps 


(өзі) А 
ху xp a(x}, x?) 


In the general case we want to solve a system of equations of the form 
f(x) = 0, (1) 


where f(x) = ( A(x), (х). ..., Лбх) and x is the N-tuple (x,, x3. ..., ху). Equation (1) 
can be expanded using the first-order Taylor approximation to-give the iterative 
scheme 


xt 1) = о _ 7". 1 f(x?) (2) 
where J is now given by 
[ Hoy... Ae 
x, Oxy (x) 
oy... Ar 
ax, Ixy (x) 
A= Я 


Sa... Tiny 


ex, Ox; Oxy 


which is a direct generalisation of the two-equation case. The determinant of the 
matrix J, is known as the Jacobian of the functions { f! at x. Throughout this discus- 
sion we have assumed that the inverse J ^! exists, which is the same thing as saying 
that the Jacobian, det[J,], is nonzero in a neighbourhood of the root x = a of Equation 
(1). Actually this ensures that we can find a neighbourhood of x — a in which « is 
the only solution. 


In using the scheme (2), an inverse need not be computed at each iteration: instead a 
linear system of order N has to be solved. We can see this by noting that Equation (2) 
can be written in the form 


J gem xU? — х'"* D) = хо) 


which is a linear system which сап be solved for the vector x) — x@* U, 


Although we shall not investigate this method any further it must be stated that the 
choice of initial guess is important with Newton's method because if the approxt- 
mation is not close to the solution, the iteration will not converge in general. Fortun- 
ately in our problem we can often get a good initial approximation by extrapolation 
from previous levels. 


Example 


Evaluate the finite-difference solution, along the first time level, of the nonlinear 
problem 


2112 
oa 0<х<І, 1>0, 
U(x, 0) = sin zx 0<х<і|, 
U(0, г) = U(1, 1) = 0 і>0, 


taking г = k/h? = 1 and h = 4 in the implicit scheme given by 
k 2 2 2 2 2 2 
Ші Wig — gp быа = 2и + Wiger t Uri 2и] + ui igh 


and solving the nonlinear equations by Newton's method. 


Solution 


The region of interest is shown in the following diagram. 


th 


2 


i=0 i=l i= i-3 i24 x 


- There are three unknowns along the first time level, и, |, u; ; and и; ү. By symmetry 
нул = из, SO we need only consider iu, | and и; , which we shall write as u, and uz 
respectively. Substituting the initial and boundary conditions into the. difference 
equation yields two equations for и, and и, given by 


fiu, u2) = —2uz u$ — 2u, + J2 = 0, 
Ди, из) = 2u? — 2u — 2u + 1 = 0. 

Using the form of Newton's method given by 
Jalu — Ш 1 = f(a), 


where u = (u,, u2), we have 


—4и—2 210) um — uot u — 2(u? + (uy)? 5% 24% EE J2 
AW tl <2] |7 (uly? — 294 — 2u +1 |, 
since 
9f, Әл д 
= -4,-2, T5224, 22 = da ens 
ди, і диҙ Ші Qu, ы Qu; ds 


Starting the iteration with (ufP, uff?) = (t; о. и» 9) = (24/2. 1) we get 


[—48284  20][0.7071 — ut" 0 
2.9284 —6.0 | | 1.0000 — иу) -2 
from which (ut, ut!) = (0.5355, 0.5857). The second iteration yields 
Г-41420 1.1714] [0.5355 — ug 0.1127 
| 2.1420 —4.3428 | | 0.5857 — ug — 0.2840 
and so (u'?), и?) = (0.5456, 0.5253). The third iteration gives 
а 1.0506 | | 0.5456 — и? 0.003438 
2.1826 — 4.1012 | | 0.5253 — ш?) — 0.007075 


from which (и), ut?) = (0.5461, 0.5258). This final iteration leaves residuals 
(—7.622 x 1075,1.379 x 1075). Hence, 


Uy) = U3, = 0.546, 
Изд = 0.524. 


SAQ9 


Evaluate the finite-difference solution along the first time level of the nonlinear 
problem 


ға - = 0<х<1, t>o 
U(x, 0) = 16x(1 — x) + 1 0<х<1, 
U(0.0) = U(1,) = 1 tz 0, 
Taking r = 45 and h = 4 in the Crank-Nicolson type formula quoted in the example 


above. Perform only the first three iterations using Newton's method. 


(Solution on p. 30.) 


15.1.4 Three-Level Schemes for Nonlinear Equations 


We have seen that finite-difference methods based on two time levels for both linear 
and nonlinear equations have poor stability properties if explicit, and are difficult 
(complicated) to use if implicit. Therefore, we turn to three-level schemes which, in 
the linear case, can combine good accuracy with good stability properties although 
they are not very economic in computation time. We shall now proceed to investigate 
the nonlinear parabolic equation given by 


wu) ee = 2 (о) with aqu) > 0, btU) > 0. () 
ôt дх дх 
The simplest three-level explicit approximation to this equation is given by 
Up jay — Wij- 1 
ДОЯ) БЕ T = тз Slalu 8... (2) 


There seems little point in pursuing this scheme further because we have already 
seen in the example in Section 15.1.1 that when a(U) — b(U) — 1 it is completely 
unstable. However, in the linear case we obtained unconditional stability by using 
the Crank-Nicolson device. This fact indicates how we can proceed here. Expanding 
Equation (2) gives 


bu; (ува — j-)) = Rr {a(i (uis ij = ш.) = абша.) — ti- a 


where we have used the central-difference formula 


ди; = Ue ye 


А VIheit-averiees oq ime levels j — 1. 
We now replace ш, у. u; and шү; by their averages over three time le і 
j and j + 1. so that 


1 is replaced by Мира + eau + наа 


VENT 
u; is replaced by (и. + tij + tija) 

and 

ij. is replaced by ліш; i jsi + cau + oui 

and we replace +; j by an average over neighbouring mesh points so that altis, j) 

becomes 


PAT [eos te) 


- 


This process leads to the finite-difference scheme given by 


and a(t, j) becomes 


зр 
bln ужа — usd irt ped тырна Шыт Ug b Ue au = Ui j- 1) 
— Wu ie. Z Mien gen c UL; T inaug js. — ша]. 


This replacement does not alter the order of accuracy which remains O(h?) + O(k?). 
Note that for nonlinear equations Lax's Theorem cannot be used and we have to 
aim directly at a proof of convergence. In fact, it has been shown that this scheme 
is both stable and convergent for all values of r. The most remarkable feature of the 
scheme is that, although implicit. it gives rise to a system of linear algebraic equations 
to be solved at each time level. Hence. the above scheme avoids the complication of 
solving sets of nonlinear equations inherent in the Crank- Nicolson method on two 
time levels! 


The following SAQ looks at the accuracy of the above scheme. The solution is quite 
involved since you will need to expand арапа 2; by Taylor's Theorem. 
SAQ 10 


Show that the local order of accuracy of the above three-level scheme when applied to 
the nonlinear equation (1) is O(h?) + O(K?). 


(Solution on p. 30.) 


PDE 15.2 


15.2 PARABOLIC EQUATIONS IN TWO SPACE DIMENSIONS 


A typical parabolic equation in two space dimensions is the two-dimensional heat 
equation 

oU QU PU 

a d ys г> 0, x re(0, 1). (1) 
We shall examine the numerical solution of this equation in a square region of the 
xy-plane in the usual way by covering the region with a rectangular mesh with 
spacing h in both the x- and y- directions. The mesh points (x;, y;.f,) are given by 
x; = th, y; = jh and t, = nk where i, j and n are integers and i =n = 0 15 the 
origin. The solution at (x;, у;,1,) of a given finite-difference scheme will be denoted 
by t; jn» and the true solution of the differential equation at that point will be denoted 
by U(x;, yj, ta) or more simply, by U 


Ыл 
The standard explicit finite-difference replacement of Equation (1) can be obtained 


by replacing the time derivative by the usual forward-difference formula and the 
two space derivatives by the usual central-difference formula. Thus, we obtain 


1 е5 

ptus = т (6; + (ЭТТ (2) 
ог 

Wi pnt = Mijn + (Шаан + rau Шал + Шс Au; junh (3) 


where r = k/h? is the usual mesh ratio for parabolic equations. 


SAQ 11 
Find the local truncation error of formula (2). 
(Solution on p. 32.) 


The stability of a formula used to approximate equations in two or more space 
dimensions can be investigated by von Neumann's method which we discussed in 
Section 8.4.3 of Unit 8, Stability. In the case of Equation (1), where there are three 
independent variables, we look for separable solutions to the error equation which 
have the form 


E 


pan = eif PhgiBsnhgn. 


The error equation is obtained by considering starting the computation with the 
vector ug instead of ug when п = 0 and putting 


E = 1 


p.g.n рал > l 


LE 

The von Neumann criterion for stability is then 
[sn 

as is easily verified by the techniques of Unit 8. 


SAQ 12 
(a) Use von Neumann's method to investigate the stability of formula (3). 


(b) Give a sufficient condition for the convergence of formula (3) when applied to 
the differential equation (1) as h, k ~ 0. 


(Solution on p. 32.) 


As usual, poor stability properties leading to increased computation mean that 
explicit methods are rarely used to solve initial-boundary value problems in two or 
more space dimensions; implicit methods are almost always used in preference. As 
we have seen repeatedly, implicit methods require sets of algebraic equations, which 
are linear if the differential equation is linear, to be solved at each new time level. 


If we now apply the Crank-Nicolson averaging technique to Equation (2) we obtain 
the finite-difference scheme 

(1 2092 — 372), рл = 0 3rd? + brô?) jn- (4) 
As we are dealing with two space dimensions, we effectively have to solve a boundary 
value problem at each time step. We can therefore exploit our knowledge of the 
numerical solution of such problems (gained from Unit 11, Boundary Value Problems) 


as the following SAQ illustrates. 


SAQ 13 

If U is specified on the boundaries of a square region0 <x < 1,0 < y < 1, fort > 0, 
show that the SOR method, discussed in Section 11.4 of Unit 11, can be used to solve 
the system of equations arising from the use of Equation (4) to approximate the solution 
of the differential equation (1). Find the optimum relaxation factor when there are 
M + 1 mesh points in each of the x- and y- directions. 


(Solution on p. 33.) 
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15.3 DISCONTINUITIES IN THE SUBSIDIARY CONDITIONS 


The very nature of finite-difference approximations to partial differential equations 
means that they are unsatisfactory near. discontinuities. We. have usually derived 
finite-difference replacements to derivatives from: first- or second-order Taylor 
approximations which we know to be of dubious value at or near discontinuities. АП 
the estimates of local order of accuracy and global error for such schemes depend on 
the boundedness of partial derivatives of the true solution to the given partial differen- 
tial equation. These estimates may not be valid-if.the solution has discontinuities at 
one or more boundary points. Obviously great care has to be exercised in these cases 
and usually we have to treat each case individually. 


The most common form of discontinuity in initial-boundary value problems occurs 
at the corners of the solution domain, where the initial line meets a boundary line. 
If the initial line is t = 0 and the other-boundaries are x = 0 and x = 1, typical 
initial and boundary conditions are given by 


U(x,0) = fi) UOD = ft). UL) = fs, 
which we usually take to mean 

im U(x, t) = fi(x) 0<х<І, 

Jim, U(x, t) = ДЫ) t>0, 

Jim U(x, t) = f(t) t>0, 


respectively. The corners (0, 0) and (0, 1) sometimes involve some form of discontinuity. 
The worst type. of discontinuity is in the function value where, for example, 


ЛО) = 500) ог AC) е f0). 


Weaker forms of discontinuity also exist. For example, suppose that we are con- 
sidering the diffusion equation, with initial and boundary conditions of the above 
type, given by 


fit) = 0 
fold) t. 


Obviously f,(0) = f,(0) = 0 and there is no discontinuity in the function values at 
the point (0, 0). However 


The effect of such discontinuities does not penetrate into the solution domain for the 
analytical solution. However, this is no longer true for finite-difference solutions 
which have doubtful value in the neighbourhood of the points of discontinuity. 
Usually we do not have to worry because the true solution becomes "smoother" as t 
increases and any errors introduced in the finite-difference solution decay, provided 
the methods used are stable. We have seen a good example of this in Section 5.2 of 
Unit 5, Initial Value Problems. However, this is not always the case, as you will see 
in the next reading passage. 


As a simple example of how analytical numerical methods can be used together to 
advantage we consider Exercise 2 оп S: page 46. There is a discontinuity at the 
origin of the given problem which means that the finite-difference solution is inaccurate 
for small values of x and t. (See the table in S: page 48.) If, however, we take the 
analytical.solution at t — 0.001 (given in Table 2.28) as the starting point for the 
simple explicit scheme using r — 0.1 (given at the top of S: page 47) we obtain the 
results in the following table. 
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in 


x=0 0.1 02 0.3 04 0. 


г= 0001 | 0 0.9747 1.0000 1.0000 1.0000 1.0000 
0002 |0 0.8798 0.9975 1.0000 1.0000 1.0000 
0.003 | 0 0.8036 0.9860 0.9997 1.0000 1.0000 
0.004 | 0 0.7415 0.9691 0.9984 1.0000 1.0000 
0005 | 0 0.6901 0.9493 0.9956 0.9998 1.0000 


If we concentrate on the results at x = 0.1 where the finite difference results are 
poorest we find that, at г = 0.002. the percentage error was previously 7.5 per cent 
whereas now it is down to 0.7 per cent. Similarly, at г = 0.005 the percentage error 
has improved from 3.8 per cent to — 1.1 per cent. 


M 
is) 
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15.4 SUMMARY 


The quest for simple schemes with good accuracy leads us to look at three-level 
schemes. We have seen how to adapt the matrix method to investigate the stability 
of these schemes. 


We have investigated a simple class of nonlinear partial differential equations and 
showed how our idea of stability has to be amended to local stability before we can 
investigate our simple finite-difference schemes analytically. We have seen that, in 
general, implicit methods for nonlinear equations necessitate the solution of sets 
of nonlinear algebraic equations and we presented a straightforward method for 
their solution. However, this method illustrates the major difficulty in solving non- 
linear algebraic equations, namely the large amount of computation which is 
necessary. Fortunately. for one simple class of problems, we were able to quote a 
three-level implicit scheme which gives rise to linear algebraic equations. 


We have also looked briefly at the diffusion equation in two dimensions. Finally, 
we have taken a look at the problems associated with discontinuities in parabolic 
equations and made the important point that finite-difference methods should not 
be used in isolation. It is important that numerical methods should be used in con- 
junction with analytical techniques in order to obviate those difficulties which usually 
arise in real problems. 


15.5 SOLUTIONS.TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 
The local truncation error of the three-level scheme is given by 


Uii; = Vigna _ Vieng 7 205 + Uii 


Tij 2k k? 
where U is the true solution of the differential equation 
au OU 
д 0x? 
Taylor's Theorem gives 
Uy jar — Шш = 2k 204 + Башы + 0(к5) 
and 
Ui 1; -2U;; + Vinny =? PU * UM + (nS). 
Hence, 
Tj = O(k?) + O(h?) 
since 
au PU 
DEP 


Thus, the three-level explicit scheme is locally more accurate than the two-level 
explicit scheme, for which, as we saw in Unit 5, 


T, = O(k) + O(h). 


Solution to SAQ 2 
The Du Fort and Frankel scheme can bé written as 
(L + 20) уу = Ariana + ш) + - 20u;-1, 
that is, 
Аша = Bu; + Cuj., 
where A = (1 + 271, C = (1 — 2r)I and B is given by 


Ot 
1 01 
1 0 1 


The square matrices A, B, C and I have order M — 1, where M — 1 is the number of 
internal mesh points along each time level. The equation 


det324 — åB – C) = 0 
now becomes 


122 
ЕЛ8 


Й B 
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and if и is an eigenvalue of B then А is an eigenvalue of the amplification matrix if it 
satisfies 


(14-27022—(1— 2r) 
Ee 

À 
The eigenvalues of B are given by (see Appendix to Unit 8) 

5л 
= 4p Bus = E - 
Hs r cos А 5=1,2,...,М – 1, 

and hence 


(1 + 2722 — 4r cos (5) 4, =- (1— 2r) = 0, 


from which 


For stability we require 


max |A*| < 1. 
s 
If 4r? sin? (5л/М) < 1 this becomes 
x 
—(1 + 2r) < 2r cos (ш % | — 4r? sin? ШІ <1-?2ғ, 
which is clearly satisfied. 


If 4r? іп (вл/М) > 1, we have r > 4 and 


е с H 
nc 2r cos іш t (е sin? іш - ] 
ai ык акан т шышаны 


142r (Ге pen 
so that 
4r? cos? А + 4r? sin? г -1 
£92 a 
el = (1 + 20 
_ 42-1 
| (1 + 2r? 
E 2-1 
142r 


and | 4 < 1 for all r > 4. 


Hence, the scheme is stable for all (positive) values of r. 
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Solution to SAQ 3 
(a) The order of accuracy of the Du Fort and Frankel scheme is 
O(h?) + O(k?) + Ol? ih): 
we can plot the order of each term against p where k = O(h”) as h ~ 0. 


order 
a 


OK?) = O(h?) 


OU?) = 008772 


Oth?) 


A 


The order of the method is shown by the dark line. Thus the maximum order 
of the method is O(h?). 


(b) We want to obtain the largest order in the method and at the same time achieve 
the largest step size k (as h ~ 0). Now since h is usually small and К = O(h’) 
we require p to be the smallest value such that the order remains O(h?). That 
is, p = 2. 


Noticethatp = 2meansthatr = k/h? = constant. We prefer to haver = constant 
in practice since this has the added advantage of giving convenient calculations. 
Solution to SAQ 4 


The local order of accuracy of the Crank-Nicolson formula is given by O(h?) + O(k?). 
Putting k = O(h?) as h ^0 gives the following graph of order plotted against p. 


order 


O(K?) = oue) 


O(h?) 
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SAQ 4-5 


The order of the method is shown by the dark line. Thus the maximum order of the 


method is 2. Since h is small, the largest value of k is obtained when р = 1. That is, 
we prefer to take К = O(h). 


Solution to SAQ 5 


The following is a sample output from the program $COM2321 which uses both the 
Crank-Nicolson and Du Fort and Frankel schemes to solve the given problem. 
The program has one input parameter N which is the total number of mesh points 
along each time level. The results from the schemes are compared at t = | by printing 
out the errors (the difference between the computed solutions and the analytical 
solution). The output consists of two tables of numbers. The first is a comparison 
between the Crank-Nicolson scheme with k = h and the Du Fort and Frankel 
scheme with k = h?. The second table shows the errors of the Crank-Nicolson and 
the Du Fort and Frankel schemes both with k — h. 


RUN 
COM321 


A PROGRAM TO COMPARE THE EFFICIENCIES OF THE CRANK-NICOLSON 
AND DU FORT-FRANKEL SCHEMES 


INPUT THE NUMBER OF MESH POINTS ALONG EACH TIME LEVEL (N)?7 


SOLUTIONS EVALUATED AT Т-1.И 
TIME STEPS) 
TIME STEPS) 


CRANK - NICOLSON SCHEME WITH K=H ( 6 
DU FORT - FRANKEL SCHEME WITH K-H*2 ( 36 


ERROR IN 
X CRANK = М DU FORT -F 
.166667 2.5й243Е-И85 2.58457Е -05 
.333333 4.33536Е -05 4.4766ЙЕ-й5 
45 5.0й0798Е-05 5,.16914Е-й5 
.666667 4.347U4E -05 4.47661E -85 
833333 2.48565Е -45 2,58458Е-И5 
DU FORT - FRANKEL SCHEME WITH K-H ( 6 TIME STEPS) 
ERROR IN 
X CRANK — М DU FORT -F 
. 166667 2.50243E -95 319499 
. 333333 4.33536E -И5 .553387 
.5 5.00798 -45 .638996 
.666667 4.34704Е-й5 .553387 
.833333 2.48565Е-И5 .319498 
DONE 


We have used a value of 7 for the input parameter N., giving h = 4 for both schemes. 
The first table shows that, as predicted in the text, the Du Fort and Frankel scheme 
with k = h? gives an accuracy comparable with that obtained using the Crank- 
Nicolson scheme with К = h. [With N larger than about 7 it will be found that there 
is a considerable waiting time between the computer printing the message 


DU FORT - FRANKEL SCHEME WITH К = H T2... 


and the table which follows it. This time is mostly taken up by the computer evaluat- 
ing the Du Fort and Frankel scheme. There is no noticeable time taken in computing 
the Crank- Nicolson solution.] 
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The second table. in which k = h for both schemes shows clearly the inferiority of the 
Du Fort and Frankel scheme: the errors do appear to be extremely large. We can 
account for this by reference to the expression for the local truncation error of the 
Du Fort and Frankel scheme, 

к? BU; ; Е h? 2*%0,; z k 020,; 
8-6 aD 12 ә '"mw д? 


where we see that the ratio k?/h? is unity when k = h. This means that the local error 
will contain large contributions from the term @?U, ,/dt? and it is likely that the com- 
puted solution is a better approximation to the hyperbolic equation 


2 2 
QU PU PU у 


à 22 * а? 
than to the given parabolic equation. (See SAQ 7 of Unit 8.) 


Solution to SAQ 6 
In full the three-level implicit scheme is 


2r 4r 2r 2r 4r 2r 
Ui jl 


=F Misije + |: + 


3 т-н = шж g gj + 34-4 


3 


2r 4r 2r 
+ 3-і * f = 5 tij-1 + 3 agnis 
This scheme gives a tridiagonal set of equations, which can be solved by the recur- 
rence form of Gauss elimination just as for the Crank-Nicolson scheme. The only 
differences between this scheme and the Crank-Nicolson are in the coefficients 
(which are again constants) and in the evaluation of the right-hand sides. Now we 
need more computation to evaluate the right-hand sides of the three-level scheme 
than in the Crank-Nicolson scheme. Therefore, on the grounds of efficiency, we would 
prefer the Crank-Nicolson scheme. 


Solution to SAQ 7 


Expanding the formula, we obtain 


—2rui 1,44 (3 + Ф)ш jay — жыра = 4u,; — Ші-і» 
which has the following molecular diagram. 
Еј i+ jet 


ij-1 


Once again a tridiagonal system of equations (with constant coefficients 
solved. This scheme requires less 


does the Crank—Nicolson scheme 


| tior ‹ ) must be 
arithmetic in evaluating the right-hand sides than 
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Solution to SAQ 8 


The Crank-Nicolson type formula for the equation 


au Хи” 
д Ax? 
is 
u -и 
ем ра = gp Uo asi — 2upaei Ша T pei — дин а + Me iq) 


The introduction of an error E,,, into и, given by 
E, 


и и 


* = 
P.4 P 


leads to the equation 


та — 


* = — (uk. — = 15" - жт 
(ира+а Up 1) (ира ира) = [uz lal Up 14817 2(ира+1 mL 
жт cuui жт 
+ up gta T 5-і + Ubri — Ир 
= жт _ т km __ 
2454 ира) ирлан up- 1.a): (1) 


Assuming that и varies slowly over a small region in comparison with the errors and 


that E, is small compared with u,,, we have 


Up+iqg+i S Upgtr ~ Чр ража US qu = Upg US 
and 


m 


"—1 
pq тЕ, pa 


жт __ 
Шра Ц pal p 


Equation (1) now becomes 
Epari = Epa + тира (Еу ду — 2Ё ду ЖЕ, рө + Ёл 7 2Epq Epig. 
Substituting the relation 

E 


= gilphzaq 
pa = ё "©, 


we have 


E=1 + ти" (ele — 26 + e^ PE + eh — 2 + gc ith) 
on dividing by еї?" Putting R = rmu7;! in this equation and rearranging, we 
obtain 

h . 3 Bh 

ё|1+ 2R sin? P" =1- 2R sin? P^, 

2 2 
and clearly |é| < 1 for all values of r provided u7; ! > 0. The method is therefore 
locally unconditionally stable. 
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Solution to SAQ 9 


тА 
А 1 1 
j=2 
Я 1 uy и; и, 1 
pel 
I 
jo 4 5 4 x 
і-0 #= | i=2 i23 і= 4 x 


Let the pivotal values on the first time level Вен, = и, = изу and u, = и. The 
nonlinear equations which have to be solved are given by 


Silty. и) Bad — 2u} — 20u, + 75 = 0 


and 


Differentiating we obtain 
6f, of, @ GE 
LA = —4u, — 20, а _ 21. h _ 4u. 22 = —4u, — 20. 
ді ди» Qu, ди, 


Starting the iteration with (uO, ut?) = (и, o, моо) = (4. 5). we obtain 
[-36 101Г0.656 _| 712 

| 16 —40]]1162] | —36 

and so (i! 1) = (3.344, 3.838). The second iteration gives 

І [-3334 767517001234] Г 0.8028 

| 1334 —35350]| 005098 | ~ | —1.9667 

2) 


from which (12), uf?!) = (3.356, 3.787). The third and final iteration is given by 
| [ —3342 7.573] [9.325 x 107? _ | 703107 
L 1342 —35.146 || 1.346 x 107* |. | 0.1205 


and hence (uf), ы) = (3.347, 3.786). 


А convenient way of tackling this problem using the computer would be to use the 
library program $GEM201 which solves the matrix equations by the Gauss Elimina- 


tion Method. The coefficients could be obtained from a simple program such as the 
following. 


10 INPUT Xi. X2 

20 PRINT X212-2«X112—20«X1-- 75: 2«X212--2«X112 —20«X2 1.82 
30 PRINT —4« ХІ – 20:2+ X2:4« Xl; —4« X2 — 20 

40 STOP 

50 END 


Solution to SAQ 10 
The local truncation error is given by 


bij j 
р x Un =4 U; 5-1) 
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1 
spi Uni сач Uj jet + О. 09; gt Ui. i- i" U;;- ) 


= Шужа — Шужа + Ui — Uii + Фуа — Ха-а] 


where b; j represents b(U; ). By Taylor's Theorem 


1 д0, ; 
zg Cui —-U,.)- x + O(R, 


1 
sii Vis itt 7 Vases + Una — U;; + Ui.ij-, = Uij- 1] 


[ата no 
~ [лдх 20x? " 60 


+ ото + O(k?)) 
and 


I 
320+: = Ора + 9;- Ui- il 9;- к Ш-аш-1] 


12 129 n? Е 
-|-—-—- k? 
ЕЕ sat ca 96 JONES ). 


Now, by repeated use of Taylor's Theorem, 


йу Uig 

өй = «[ ыы 5 "ы 
һд0,, Һ?д?Ш,, 
к. ax + % Ox? 
4 Utoy ÈU day dU)? ёа, 
T9 ôx dU 4 0x dU? 


+ “Л 
+ O(h?) 
апа 


2 
h0U;,da,, WP PU sda; | I? y te 


= Өф. э ае aj al. TUB h3), 
Se dU * 4 oe dU (х) dui ^ QUÀ 


= md Е Ж des) 


where а; j represents a(U; j). Hence 


P 1820, 
Ty = by al + Ой?) — (o£; + ор) [5 air o0 + ow 


тр ES + об?) + ow] 


‘J 4 O(h?) + O(k?). 


ди) PU, e J da; 
j dU 


QU ð 20 
60) == TO («nz 
да QU ди 
Saa | up 
ди)? da д? 
= | 19 + alU) т" 


and so 
Tj = O(h?) + O(K?). 
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Solution to SAQ 11 


The local truncation error T; ; , is given by 


быша — был 
Tja = мы ты 


Ui iat Ui-i ja t Viggen t ÜUij-in 7 4U; js 
Taylor's Theorem gives 


h? 
ду, KOU jn | k? Ui ja 
Lp. ijan ji ders 
Ui; ei 7 Uij =k er +5 EE 6 аз 
ди,” Ui i, 
Uisi ja t Ui-ija = ijn + PA 12 axt 
and 
QU, В 9*U,i, 522 
Орла + Uij-in = 20, +h? or a Do + 
Therefore, 
Win PU in OU, + кё in 
Mate Ax? oy ^2 08 


h? д*0, jn + 9*U,;, digi 

12| Ox* , 

The differential equation gives 
ди ди ди =0 
й Ox! ду | 

and so 


Тұл = O(k) + O(h?). 
We see that 


lim Ty jn = Qash, k ~=0 


and the given formula is consistent with the two-dimensional heat equation. 
Solution to SAQ 12 


E 


we perform the calculations with u* replacing u. and 
p.q.n ши 


(a) Suppose that we introduce an error at each point along a time level such that 
pan T u 


* 

pnt 

The error E, then satisfies the difference equation given by 
Putting Ep qn = efiPheilishe we obtain 


Еа = Ё „+ (Ej, tan ЖЕ, + Ера Epa-1a — 4E, 


pagn) 
€ = lt r(eP 4 ейән 4 -ibih | eripi _ 4) 


after dividing through by еі'"е! 2" Since 


ет + emiBh 2 cos Bh 


25 agen? {Bh 
— 2 — 4sin e | , 
we obtain 


cce etl 


Now, for stability || < 1 and so we require 


-1«1- ur (2) + sin? ЕЛ <i. 


PDE 15 SAQ 12-13 


The right-hand inequality is trivially satisfied for r > 0 whilst the left-hand 
inequality gives 


2M Bah sin? z 


and since 0< < sin? $f,h < 1 and 0 < sin? 3f;h < 1 we get the stability condition 
"< ў. 


rs 


(b) As the scheme is stable if r < 4, and is consistent with the given differential 


equation, it is also convergent by Lax’s Theorem. (Note that both the differential 
and difference equations are linear.) 


Solution to SAQ 13 


The finite-difference method gives 
Mijnel — 24» LjasiCt Mi-ijaer FM jetta + Ujj-iae] — 4и jinsi) 


— known quantities, 


ie. 


А 
Шіл ~ 20725 120) (калжа F Menge + Шужілғі H Ujin) 


= known quantities. 

The collection of all equations for 1 < ij < M — 1 can be written as 
Au =b. (1) 
This formula has precisely the structure of the five-point Laplace formula discussed 
in Unit 11; therefore the discussion about SOR applies provided that the associated 


Jacobi iteration converges. The natural ordering (Section 11.4.2 of Unit 11) leads to 
the following form for the matrix A in Equation (1): 


B 
I 


where A is square of order (M — 1)?,@ = —r/[2(1 + 2r)], B isa square matrix of order 
M — 1 given by 


011 
1- 9-51 
100781 


I 2021 
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and Г is the identity matrix of order M — 1. The associated Jacobi iteration matrix J 
is given by 


D I 
I D 1 
1 DI 
Ј=1-А= –0 
I D 1 
I D 
where 
[o 1 
1 0 1 
L0 3 
Dz 
1.0 1! 
1 0 
is of order M — 1. By the result of the appendix to Unit //, the eigenvalues of J are 
given by 
r л л 
фат T+ ze А + cos (27) ] ьа = 1,2,...,М – 1. (2) 
Since 


рп qu 
1 «cosy, cos sr < 1 Pq = 1,2,...,М – 1 


we see that 
—2r 2r 


1:2: 5%а<Тұ% 


and therefore |a| < 1 and the Jacobi iteration is convergent. Тһе SOR method is 
therefore convergent too. To find the optimum relaxation factor Wap We use the 
formula 
ie 2 
om T+ tl = py} 


as given in Section 11.4.3 of Unit 11, where p(J) is the spectral radius of the Jacobi 
iteration matrix. Now, from Equation (2), 


cos = 
M 


2r 
00) = TT, 


giving 
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16.0 INTRODUCTION 


It is unnecessary to stress the importance of understanding the flow of blood in the 
circulatory system. This system is extremely complex and it is unlikely that a satis- 
factory mathematical theory of the whole of it will ever be obtained; indeed, at 


present we are some way from obtaining an accurate mathematical model of arterial 
flow. 


This unit describes some of the work done with a view to understanding the formation 
of thick plaques in the arterial wall which partially block the artery. This disease, 
called atherosclerosis, is the chief cause of male deaths in Western societies (New 
Scientist, 12 July 1973, p. 65). The reason for the formation of these plaques is not 
understood, but sometimes the artery becomes sufficiently restricted to warrant 
remedial surgery, which may consist of the insertion of a by-pass; this is either an 
artificial artery or a segment of the patient's vein. Unfortunately it is sometimes found 
that both the parent artery and the graft tend to block again at the down stream end. 


There is evidence that, in the formation of plaques, the chemical composition of the 
blood is an important factor, but now physical properties of the flow are also con- 
Sidered important. At present it seems that a combination of chemical and physical 
effects are important in the formation of plaques. The work described here is an 
attempt at understanding some of the physical properties of blood flow. 


In order to obtain some understanding of flow in an artery there are three possible 
approaches. Firstly, experiments can be performed on animals; secondly, we can make 
models in the laboratory; finally, we can make mathematical models and experiment 
with them. In practice the first option is not entirely satisfactory for detailed observa- 
tions, because of the small diameter of the arteries in small animals and because of the 
difficulty of doing accurate experiments in situ, although such experiments are 
necessary before the other two methods can be used* The mathematical models are 
usually so idealized that, without further evidence, their results cannot be relied upon 
to give accurate quantitative results. One of the problems is the very difficult task of 
gauging accurately the errors produced by making the approximations. However, 
the theory can give an insight into qualitative relations and into the likely important 
parameters. Physical models can be made of systems too complex to be modelled 
mathematically, and can give some idea of the significance of the approximations 
made in the mathematical model; but these models are subject to their own errors. 
In practice a combination of all three methods is necessary to make progress. 


In this unit we shall consider an elementary mathematical model, but where possible 
we shall compare our results with experiments. Before considering blood flow in 
more detail we consider some general properties of viscous flow (Section 16.1), and 
in Section 16.2 we shall derive some equations pertaining to flow in cylindrical tubes. 
In Section 16.3 we consider in more detail the approximations made; and in Section 
16.4 we apply the work of the preceding sections to blood flow. In Section 16.5 we 
compare our results with experiments. 


Note that in this text we do not restrict ourselves to the customary SI units; this is 
because most of the literature on this subject uses the cgs system. The symbol д is 
used for two distinct purposes; the first is to denote the coefficient of viscosity of a 
fluid, and the second is to denote one micron (= 107 $ m). There should be no confusion 
between these two uses, both of which follow convention. Finally, following the usual 
practice, pressures are measured in millimetres of mercury (mm Hg); since the density 
of mercury is 13.6 gm/cc, the pressure needed to support 1 mm Hg is 133.4 Nm 
or 0.019 Ib per sq. in. 


* A rat will be mailed in the next package. 


161 LAMINAR FLOW IN VISCOUS FLUIDS 


Throughout this unit we shall be concerned only with laminar flow, that is, flow in 
which the fluid moves as if it were in layers with different velocities. Not all flow is 
laminar, as can be seen by observing the passage of an object through water; if the 
velocity is small the fluid flow will be laminar, but if it is large enough eddies will 
form and the motion will be turbulent. In the case of turbulent flow the equation 
of motion is nonlinear. so a mathematical description of the flow is very difficult. 
In general, however. for the case that we are considering in this unit the flow is laminar. 


All liquids are viscous, but in some cases the viscosity may be neglected: whether it 
can be neglected depends not just upon the liquid but also on the circumstances 
of the flow. In blood flow conditions are such that viscosity cannot be ignored. 


If portions of a liquid are caused to move relatively to one another, the differential 
motion gradually subsides unless sustained by external forces; conversely, if a portion 
of a mass of liquid is kept moving, the motion gradually communicates itself to the 
rest of the liquid. We require a mathematical description of these effects, which are 
cvident from immediate observation. 


x 


Newton was the first to formulate a hypothesis regarding viscous forces. His hypothesis 
amounts to the following: if two parallel plane laminae of area A move in the same 
direction with different constant velocities v, and р,, parallel to their planes, the force 
required to maintain this differential motion is. 

HA 

T (01 — 05), 
where Az is the (perpendicular) distance between the laminae. Since the velocity of the 
fluid changes continuously we may replace the differences by a derivative, so that 


— 0 


The constant y is called the coefficient of viscosity; it has the dimensions ML-'T-! 
and in the cgs system the unit of viscosity is called a poise, in honour of Poiseuille. 


In practice, the convenient unit is one hundredth of a poise or a centipoise (cp); the 
viscosity of water is about 1 cp. ` 


Newton did not pursue the subject very far (although he did consider the case of 
flow between two concentric cylinders), nor was the problem studied again for more 


PDE 16.1 


than a century. Nevertheless this first contribution is commemorated in the use of 
the term Newtonian fluid for simple viscous fluids. A Newtonian fluid is one whose 
viscosity does not depend on the velocity gradient. 


The eighteenth century produced many great mathematicians, of whom Leonard 
Euler and Daniel Bernoulli devoted much thought to hydrodynamics, but their 
work was confined to ideal fluids, that is, fluids without viscosity. At the end of the 
century Coulomb studied the damping of the oscillations of a disc suspended in 
liquids of different viscosities and made the important observation that the smooth- 
2. roughness of the surface of the disc did not greatly influence the drag of the 
iquid. 
The first work on flow in cylindrical tubes appears to be that of Girard in 1813 using 
brass tubes of 2 to 3 mm in diameter. He obtained the relationship 
PR? 

0%-1- 
where Q is the volume flow per unit time, R is the radius of the tube and P the pressure 
drop along the length L of the tube. Thus he observed that the flow varied diréctly 
with the pressure and inversely with the length but thought that it varied with the 
cube of the radius. Ten years later Navier derived theoretical equations for flow 
of viscous liquids in cylindrical tubes but obtained the incorrect result. apparently 
confirmed by Girard's experimental results, that the flow was proportional to the 
cube of the radius. 


The first published experimental work indicating that the flow is proportional to 
the fourth power of the radius was that of Hagen in 1839. He used brass tubes ofa size 
similar to those of Girard and his results were not very accurate. The exponent of the 
radius derived from his results was actually 4.12 and he assumed that the real value 
must be 4.0. Poiseuille published his first results in 1842, and a detailed paper in 
1846. On account of this paper, laminar flow in tubes is often referred to as Poiseuille 
flow. 

Poiseuille spent much time investigating the hydrodynamics of the capillary circula- 
tion. He may be regarded as fortunate in two respects. In the first place, although he 
had a training in physics, he was also a physician who wanted to apply the results of 
his investigation to an understanding of the blood circulation. Hence he worked with 
glass tubes of capillary size (between 0.14 and 0.03 mm) whereas his predecessors had 
been engineers who worked with much larger pipes. The use of smaller tubes enables 
laminar flow to be maintained more easily and also accurate measurement is facilitated 
greatly. In the second place he was deflected from his originai intention of using blood 
as a test liquid because no satisfactory way of rendering it incoagulable was known 
and he was compelled to confine his investigation to water. Blood flowing in capillaries 
shows anomalous viscous properties that would have introduced complications 
in these pioneer studies. 


Poiseuille's results are expressed by the formula 


KPR* 
РЕТ 
where О, Р, R and Lare defined above. The value of the constant К was determined 
under various conditions and shown to fall with increasing temperature. This constant 


is clearly a measure of the viscosity, but by purely experimental work it is not possible 
to relate it exactly to the coefficient of viscosity defined earlier. 


0) 


А theoretical derivation of Equation (1) may be obtained from the Navier-Stokes 
equation. The early work of Navier on the equation of motion for viscous fluids 
was amplified and corrected by Stokes in the 1840s; the Navier-Stokes equation is 
the mathematical formulation of the problem. Stokes, however, did not consider 
the particular case of flow in a tube; the solution to this problem was obtained 
independently by Wiedemann in 1856 and Hagenbach in 1860 who showed that 


P, — Pj) aR 


( 
0- L 8и 


where р is the viscosity and P, and Р, the pressures at the ends of a tube of length L. 
That this relatively trivial solution to the Navier-Stokes equation took so long to 
find is probably a reflection on the fact that few people were interested in the problem. 


We now consider viscosity in more detail. Viscosity is the outcome of two distinct 
causes: the motion of the molecules, and the action of the intermolecular forces. 
The former is operative mainly in gases, where the molecules are free to traverse 
relatively long distances (of the order of 1077 m at ordinary pressures, this being 
very much larger than the radius of a molecule, typically 1079 m) between collisions 
with one another. Suppose that we have a fluid in which two adjacent layers are 
moving at different velocities; the average velocities of the molecules in these two 
layers will differ by the same amount. This means that the molecules arriving in the 
faster moving layer from the other will have, on the average, a smaller velocity than 
those already present, and conversely, the molecules arriving in the more slowly 
moving layer from the other will have, on the average, a greater velocity than those 
already present. This tendency to remove inequalities in the macroscopic velocity 
distribution is the salient characteristic of viscosity; in the absence of externally 
imposed stresses to maintain the flow, it would ultimately bring the entire fluid to 
a state of uniform motion or rest. 


In liquids, however. the molecular mechanism which tends to dissipate the motion 
of the fluid is essentially different from that already described as operative in gases. 
Being surrounded by other molecules, the individual molecule is not able to travel 
freely from one layer to another, and actually follows a very erratic path which is not 
well adapted to the transport of momentum. There is, however, another very much 
more effective method of exchange of momentum, which depends on the continuous 
action of the intermolecular forces. Perhaps the easiest method of visualizing the 
nature of the process involved is to imagine that the molecules are connected by 
elastic strings which simulate the action of the intermolecular forces. It can then be 
seen that, if two adjacent layers in the fluid move with different velocities, each will 
tend to drag the other in such a manner as to dissipate the relative motion, in the 
absence of sustaining external stresses. This method of exchange of momentum is 
most effective when the molecules are close together, which explains why it is rather 
unimportant in gases, yet the dominant process in liquids. An immediate observation 
is that viscous forces between the layers of fluid are proportional to the area of contact. 


x 


Consider a fluid as shown in the figure above such that the fluid is uniform and such 
that each plane parallel to the xy-plane is moving at a constant velocity v(z) in the 
y-direction. If the velocity is independent of z the fluid is moving bodily, as though 
solid, and from what we have said there will be no viscous forces, except at the 
boundaries, which we shall ignore for the present. Now suppose that v(z) varies with 
z; we seek the magnitude of the force F on the plane BCDE due to the relative motion 


of the plane B'C'D'E', each of which has area A. As we have seen, the force is propor- 
tional to the area: so we consider the force per unit area, F/A. We expect this to 
depend upon the velocity gradient 00/02, so that 


-2f0- 


o , cvy f"(0) 
=| PF в 2% 

where the last expression follows from Taylor's Theorem. Since F = 0 when éu/éz = 0, 
we have /(0) = 0. Thus, to first order in the velocity gradient, 


A "à 3) 


which is the result that Newton assumed. This expression gives the force per unit 
area on a plane surface parallel to the fluid flow. 


This method of establishing a law about viscous forces may appear artificial, but in 
fact is typical of methods used for making mathematical models. First we isolate. the 
significant variables—in our case we used a physical argument to suggest that viscous 
forces depend upon velocity gradients; this is the hardest part of the exercise and 
usually requires considerable insight and experimentation. Next a law relating these 
variables is postulated—in this instance the linear relation given by Equation (3). 
The law can be tested experimentally, or it can be derived from a more general theory. 
In the former case the arbitrary constants introduced, here р, will be determined as a 
number; in the latter case they will be determined in terms of more fundamental 
constants. Viscosity is an interesting example—for gases the constant и can be 
derived from the Kinetic Theory of gases, but for liquids no satisfactory basic theory 
exists. Experimentally, the assumption that viscous forces depend linearly on the 
velocity gradient is often good even for quite large velocity gradients. However, blood 
is one of many fluids for which this is not sufficiently accurate, as we shall see in 
Section 16.3. 


In order to determine the motion of a fluid element we require the force per unit 
volume acting on an element of the fluid. 


Consider a fluid which is moving in the y-direction with a velocity which depends 
upon z, as in the previous figure. We suppose that the boundaries of the fluid are so 
far away from the region of interest that they may be ignored. We split the fluid into 
а series of horizontal layers each of height Az, so that each moves as a rigid body with 
slipping occurring only at its boundaries, rather like a pack of cards. We suppose that 
the velocity of each layer is the velocity of the actual fluid at its centre. as in the 
following figure which shows a cross-section of the motion. 


——bk ц: + Az) 


Y 


The velocity differences across the boundary planes CD and C'D' are then 


Avcp = t(z) — v(z — Az) 
Yep = e(z) — v (4) 
Arcey = v(z + Az) — (т) 


and the velocity gradient is Av/Az. so that the viscous forces in the y-direction acting 
on a slab, of area А, between CD and C'D' are 


Aten 
Fon = -uA —— 
ср Шш Де 
апа 
Aven 
Fep = pA -> 
CD i Az 


The force per unit volume in the y-direction acting on the layer between CD and 


C'D' is then 
= Fon + Fen: = Atew — Ateo, (5) 


АА: Аз Az 


Using Taylor’s Theorem in Equations (4) and inserting the result into (5) we obtain 


so that, in the limit as Az ~ 0, the viscous force per unit volume іп the y-direction is 


2% 

Е = Lx: (6) 
We must now consider the boundary conditions associated with the motion. There are 
always forces of molecular attraction between a viscous fluid and the surface of a 
solid body, and these forces have the result that the layer of fluid immediately adjacent 
to the surface adheres to it, and is brought to rest. Accordingly, the boundary conditions 
for the motion of a viscous fluid require that the fluid velocity should vanish at 
fixed solid surfaces. In the general case of a moving surface. the velocity must be equal 
to the velocity of the surface. 


The derivation of the Navier-Stokes equation for a general fluid is beyond the scope 
of this course. However, we are interested primarily in incompressible flow along 
cylindrical tubes, in which case the equations are relatively simple to derive. 


Newton's Second Law of Motion gives us (see Appendix to Unit 14, Bessel Functions), 
for a fluid of density p, 


ot 


where the operator у · grad acts on each coordinate of v and is given in Cartesian 
coordinates by 


p (s + (v ай = —grad p + Fo, 


v. grad B +, 2 + ш 
И =й E 
в "Ox — "Oy 50: 
Fo comprises the body forces per unit volume (e.g. gravity, electromagnetic forces, 
and now the viscous forces), and p denotes the pressure. 


5401 


Consider fluid flow between two fixed parallel planes, as in the preceding text, where 
the velocity v is in the y-direction and depends on z and t only. Show that v satisfies 
the equation 

Qv др m 2% +F 

LÀ up pm 

Pa фу Mat 

where F denotes the external force acting in the y-direction, p denotes the density 
and p denotes the pressure. 


(Solution on p. 29.) 
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SAQ 2 


(a) 


(b) 


A viscous fluid at rest is bounded by the yz-plane and fills the region x > 0, 
with constant pressure and no external forces. If the boundary is now moved 


in the y-direction with given velocity u(t) show that the equation of motion is 
the diffusion equation 


àv 2% 
ay їй) = "xc (9 t) x 0, 


where v = p/p, and that »(0, t) = u(t). 


If u(t) = A cos wt, determine the velocity of the fluid, and find the force per 
unit area necessary to keep the plane moving in this manner. 


(Solution on p. 29.) 


The quantity v = p/p, which will appear frequently in this unit, is called the kinematic 
viscosity. 


11 
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16.2 VISCOUS FLOW IN CYLINDRICAL TUBES 


16.2.0 Introduction 


Since we are interested in flow along tubes, we devote this section to deriving the 
equation of motion of a fluid moving in a cylindrical tube in two special cases. First, 
in order to illustrate the derivation of the equation of motion, we consider the motion 
of a fluid rotating about the axis of the cylinder. Next. we look at the more relevant 
problem of fluid moving parallel to the axis of the cylinder, and obtain a result we 
shall require in our investigation of blood flow. 


16.2.1 Rotational Flow About an Axis 


We consider radial flow in a cylinder, i.e. we suppose that the velocity is in the 
tangential direction only, and that all variables depend only on the radius r and the 
time t. To obtain the forces acting on a fluid element, we split the fluid into a series of 
concentric rings of equal width 2A, each supposed rigid and rotating with angular 
velocity equal to the actual angular velocity of the fluid at the middle radius of the ring. 
(The angular velocity is the rate of rotation, so that, if w denotes the angular velocity 
of the ring, the tangential velocity of a particle of fluid at radius r is given by ғо.) 


NY 


Let us consider the ring of inner radius r — Ar and outer radius ғ + Ar, as shown in 
the figure above. To determine the viscous forces acting on this ring it is necessary 
to find the rate at which the two adjacent rings are sliding over it. Let v denote the 
velocity of the fluid; then the angular velocity of the ring is given by 


1 


a(r, t) = 


vtr, t). 
А 


and the angular velocity of the outer adjacent ring is given by 


EC + 2Ar,1) 
"+ 2Ar,t) = -— 
on re) r+ 2Ar 


By Taylor’s Theorem, 


a2 


a (r,t) + O(Ar? 


{к 28r 1) = vtr t)  2Ar {гг} + (Ar? 
cr 


д 


апа 


2 4 
-1- 7A + gA + O(Ar’. 


1 
wlr + 2Аг, 1) = те. 1) + B(r.1)2Ar  A(r HAr)? + O(Any], 


where 


20 2r 

LEE QUE. 

rer г? 
e v 
Bez 
cr T 


The angular velocity difference at the boundary is now 
2Ar 
Aw = or + 2Ar, t) — w(r. t) = -— (Bt, t) + AG АР + O(Ar)?] 


and the relative velocity at the boundary is given by 

Av = (г + Ar)^o. 
The velocity gradient at the outer surface is then Av/(2Ar). We consider the viscous 
force acting between a segment of the ring shown in the previous figure and the 


adjacent outer ring. Let the angle subtended by the segment at the centre of the ring 
be A0: then the viscous force is . 


area x д x velocity gradient = Mr + алашы 2°. 
20r 


where / is the length of the cylinder. 


Viscous forces act on the ring, in the tangential direction, at all points around its 
circumference. Because of this, these forces may not be added directly: it is necessary 
to add their turning effects about the axis of the cylinder. Formally, we define the 
moment about P of a force F acting through Q as r x F, where г = PQ.* In our 
case we see that the magnitude of the moment of the force on the segment is 


Av 
Ur PAO —> 
(г + Ary AOp ЗА 


so that the viscous forces on the outer face of the cylindrical shell exert a total moment 
2nyl(r + Ar)? 


М; = : 


[Btr. г) + A(r. ОАР) (1) 


about the axis. 


The moment of the viscous forces due to the motion of the inner ring may be obtained 
in the same way, but it is easier merely to change the sign of Ar in Equation (1). being 
careful to note that, on this surface. the moment will be acting in the direction opposite 
to that of M,. The moment is therefore 


2nyl(r = Arn? 


М, = 


2 (Bir. t) — Atr. Ar]. 
Since the volume of the ring is 4zr/Ar. the moment per unit volume is given by 


M+M, ut {a + An? = (r — Arp 


= 53 —— "+ Mtr An? etr nnl 
4nrlAr 2r* Ar [ ] 


ap 2 


qu E + 1 To 5| neglecting terms of order (Аг). 
rer Ғ 


But the angular momentum of the ring is rrp x volume. where p is the density of the 
fluid. so that the rate of change of angular momentum per unit volume is rpér/Ct. 
and this is equal to the moment of the forces acting per unit volumet. Thus we have 


кн jj (2) 
атр 


ж For a discussion of moments. see Unit MST 282 7, Some Basic Tools. 


TI d = тг xr is the angular momentum, і = mr x f + mr x f =r x (nf) since f x = 0. Using 
Newton's Second Law of Motion. F = më. n follows that | = r x Е = moment of the forces. 
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(а) Determine the velocity of a fluid between two concentric cylinders of radii 
a. b (b > a) when the angular velocity of the outer cylinder is Q(t) and the inner 
cylinder is maintained at rest, in the two cases: 


(i) Q(t) = constant, 
(ii) Q(t) = А cos ал. where w is a constant. 
NOTE: In solving (ii) you should obtain the equation 
PV dV 
25 Ж:-- + (22 —m)V=0 
dz dz 


with m = 1. This is Bessel's equation, which has a solution J,,(z) introduced in 
Unit 14. However, being a second order differential equation, there are two 
linearly independent solutions, and another solution is denoted by Ү, (2). This 
solution was not discussed in Unit /4, since there we were looking for a solution 
which was bounded at the origin, whereas Y, is unbounded at the origin. In this 
SAQ the behaviour at the origin is irrelevant and both solutions are needed; 
you will not require any properties of Ү,. 


(b) Determine the moment per unit length on the inner cylinder in case (i). 


(c) Show that the kinetic energy of the fluid per unit length of the cylinder is 
b 
E= np [ rv? dr 


where v is the tangential velocity of the fluid and p is its density. Show also that 
the rate of change of the kinetic energy is 


* дь 
— = ap | ra, ar. 
(d) Using the equation of motion, Equation (2), and integrating by parts, show that 

dE àv "p qo? 2 

u^ mz] - 223] | (2) t ze. 


Interpret each of these terms physically. 


(Solution on p. 30.) 
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16.2.2 Flow Parallel to an Axis 


We now consider the case of fluid flowing parallel to the axis of a cylinder. 


1 


4“--------. га Аг 


т— Аг 


Consider a tube with axis along the x-axis, and in this tube a cylindrical shell of length 
1, inner radius г — Ar and outer radius г + Ar, centred upon the axis of the tube. 
Suppose that the fluid is flowing parallel to the axis and that its velocity varies only 
in the radial direction: let this velocity be w(r, г). 


SAQ4 


Show that the viscous forces on the inner and outer surfaces of the shell are 


2лші | ud 
Or orsa 


in the positive x-direction. Hence show that the total force on this fluid element is 
дер 

—4nrAr[p(x + /) — р(х)) + 2лші "3r ` 

jr- âr 


where p denotes the pressure, which may vary with time (although this dependence 
is not explicitly noted), and the shell lies between x and x + l. 


(Solution on p. 32.) 
Since the motion is parallel to the axis and the fluid is incompressible, we can use 
Newton’s Second Law of Motion directly. The mass of the cylinder is 4zrlgAr, where 
p is the (constant) density of the fluid, and its acceleration is dw/ér. 
SAQ 5 
Write down the equation of motion for the cylinder, and, by taking the limits as 
l ~ 0 and Ar ~ 0, derive the equation 
де lop и [д2 А 1де}. 
ді рдх p\ð? гдг 


(3) 


(Solution on p. 32.) 


A consequence of Equation (3) is that, if the pressure gradient is zero (constant pres- 
sure), the only possible solution in the steady state (all variables independent of time) 
is w(r, t) = 0, since the velocity is zero at the walls. Thus in order to obtain a steady 
flow the pressure must drop in the direction of flow, a balance between the pressure 
gradient and the viscous drag on the walls being achieved. We shall suppose that 
Equation (3) describes blood flow, and in Section 16.4 it will be solved with the 
appropriate pressure gradient. In the next section (16.3) we shall consider, in slightly 
more detail, the assumptions that are made and how they are likely to affect our 
results. However, before continuing, we ask you to derive a result which we shall 
require in Section 16.4. 


5406 


Determine the velocity distribution for steady flow along а cylinder of radius R with 
a constant pressure gradient — Ap. 


(Solution on p.32.) 
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163 APPLICABILITY OF THE EQUATION OF MOTION TO 
BLOOD FLOW 


Before going further and applying the equation of motion to blood flow. we shall 
consider whether our assumptions are valid in this case, and how we can test whether 
the solutions to our equations predict accurately what really happens. First we 
consider the assumptions made. 


We are assuming that blood is incompressible, homogeneous, and that its viscosity 
does not vary with the velocity gradient. The first of these assumptions scems 
reasonable and we shall not comment on it any further; the second and third are 
rather dubious. The arterial wall is approximated by a straight, rigid tube, long enough 
for the disturbance of the fluid entering it to be neglected: all of these assumptions 
could lead to serious discrepancies. 


Let us now consider blood in more detail. There are two main constituents of blood: 
there is the plasma, which is a transparent, very pale yellow liquid with a density 
close to that of water (1.03 gm/cc), and there are the red blood cells, whose major 
function is to transport oxygen. The red cells are roughly disc-shaped with a diameter 
and thickness of about 8 and 2p respectively; they occupy, on average, 40-45% 
by volume of blood. The radius of a cell is roughly that of a small capillary; con- 
sequently our fluid model could not work for capillary flow. Since the radius of an 
artery is typically 107? m (10* д) we should expect that there are circumstances for 
which the red cells may be neglected; it is known, however, that under certain circum- 
stances they can affect the flow significantly although, in general, their effect on the 
flow is not understood. 


One effect of the red cells is on the viscosity, which varies with the velocity gradient; 
the exact dependence is not known and experimental results often differ by a factor 
of 5 or more. However, in the following figure we show the general behaviour of the 
viscosity with the velocity gradient. 


Viscosity curve for human blood 


10? 

gv 

т 

а. 

9 

'* 10 
1 m ——L —— — P: a: al 
0.01 0.1 10 10.0 10? 102 10: 


velocity gradient (s ^!) 


To obtain a full understanding of blood flow it is obviously necessary to take into 
account the motion of the-red cells. This is difficult and so far has not been done. 
However, in a first approximation it is convenient to ignore them and to proceed 
assuming that the viscosity is constant. There is some experimental evidence suggest- 
ing that, in the situation that we consider, this is not such a bad approximation. 


Neglect of the red cells and their complicated effects implies that a pressure wave 
travelling in the fluid moves without change of shape; this is a consequence of the 
equations of motion being linear. But, we should expect the red cells to interact with a 
pressure wave by absorption and reflection, thereby changing its shape; the amount 
by which its shape is changed will be some indication of whether or not our assump- 
tion is justified. Of course, this distortion will depend upon the frequency of the wave: 
at low frequencies, such that the wave lengths are large in comparison with the cell 
Size, the cells will not be distorted by the wave since the pressure will be the same all 
over a cell, which will therefore not interfere seriously with the wave. At high fre- 
quencies, however, when the wave length is approximately the size of a cell, different 
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parts of the cell experience different forces, so that the pressure wave will both move 
and deform the cell, and this motion will interfere with the wave. 


To test this idea, an experiment measuring 
comparable with th: 
wave velocity 
rather larger 


the shape of a pressure wave at frequencies 
at of a heart beat {about 1 Cycle/sec) has been performed; the 
was about 1500 m/sec Biving a wave length of about 1500 m, which is 
than the size of a cell. At these frequencies no measurable effect was 
found, so that we can conclude tentatively that our assumption is justified. 


Now we consider the artery. In reality this is a flexible tube with a structure sufficiently 
complicated to ensure that its elastic Properties are nonlinear; this was observed as 
early as 1808. However, we shall have to ignore these problems; a mathematical model 
has been developed assuming that the elasticity is linear and that the tube is straight, 
and this is the next stage following the work described in this unit. It transpires 
that the effects of elasticity are small, but not small enough to be ignored completely; 


the theory presented here should be regarded as the necessary first stage of a more 
complicated theory. 


In our mathematical analysis we assume the tube to be long and to have no branches. 
These are poor assumptions; branches in particular havea significant effect, and indeed, 
as mentioned in the Introduction, the behaviour there is of prime interest. Unfor- 
tunately, at present only experiments can provide us with the information in these 
circumstances, although computer solutions to the equations of motion will doubtless 
be useful in the future. 


Clearly the assumptions that we make are severe, and can only be justified by a 
combination of physical insight and experimentation. We now consider how to 
proceed experimentally. Experiments on arteries in situ are extremely difficult and 
prone to large errors; accurate experiments are possible only if we can make a 
physical model of the artery and have all the significant parameters under our control. 
If possible we should like the physical model to be large enough to facilitate accurate 
measurements, and we consider now how the system can be scaled in size.* 


SAQ7 
Suppose that the equation 


ôw lp pfw law 

— — + -ізт +-— 
д рдх p\ôr r Or 
models the real system, that ғ and г denote the radius and time, and that the density р 
and kinematic viscosity v ( =д/р) are given. We wish to make a model larger than 
the real system and are free to scale the radius and time independently, i.e. if апа ї 
denote the radius and time in the model, we are free to choose two independent 
parameters х and f such that r = 27, t = fit. 


Show that for the model to be similar to the real system, i.e. for the model to produce 
results which can be compared with the real situation, the kinematic viscosity of the 
fluid in the model must be given by ? = 31/22, and that if the pressure gradient is 
periodic, with period T, the model period is, T/fj and the maxima and minima of the 
pressure gradient of the model are f times those of the real system. 


(Solution on p. 33.) 
For the experiment with which we compare our results in Section 16.5, the important 


parameters are given in the following table; they correspond to the choice of 2 = 0.118, 
B = 0.13. 


Апегу Моде! 
radius 0.3 cm 2.54 cm 
viscosity 28 cp 261 ср 
period 0.84 sec 6.4 sec 
pressure variation 90-120 mm Hg 11.7-15.6 mm Hg 


* Change of scale was discussed in Unit 6. Fourier Series. 
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16.4 A TIME-DEPENDENT PRESSURE GRADIENT 


16.4.1 Analysis into Fourier Components 


The solution to the equation of motion (Equation (3) in Section 16.2.2) clearly depends 
upon the pressure gradient ĉp/ĉx. and in the case when ép/éx is constant we have already 
solved the equation (SAQ 6). However. it is obvious that owing to the pumping 
action of the heart the pressure gradient in an artery is not constant. In this section we 
consider the form of the pressure gradient both in general and in the particular case 
of an artery of a dog, the latter in order that we can compare the theoretical results 


with observations. 


Pulse pressure (p) mm Hg 


Flow (Q) ml /sec 


270 


degrees of arc 


A flow velocity pulse and the arterial pressure pulse recorded simultaneously in an artery of a dog. 


The abscissa representing time in this and all subsequent curves is marked with fractions of the 
cycle length in degrees of arc. The pulse frequency was 2.75 cycles/sec (period 0.36 sec) in this 
experiment. With pulse frequencies of this order опе degree of arc represents 1 millisecond. 


The figure above shows how the pressure p and the rate of flow Q, at a point in an 
artery, vary with time through one period of the heart beat of a dog. However, it is 
not the pressure that we require, but the pressure gradient, which is determined, 
in practice, by recording the pressure at two points along the artery and subtracting 
the pressure of the downstream point from that at the upstream point at each time, 
À record of such a measurement is shown in the next figure. 


Simultancous recordings of a pressure pulse and the pressure gradient 
measured between two points 5 cm apart, in the femoral artery of a dog. 
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We shall suppose that the pressure gradient at x — 
ence at each heart-beat. Then the pressure gradi 
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Хо, Say, has the same time-depend- 
ent will be a periodic function of 


time and may be expressed as a Fourier series, so 


I). 


where c, = p, since the 
we can decompose the pres: 
of motion is linear, each h 


that we may write 
др 


дх (1) 


pressure is real, and w = 2л/Т where T is the period. If 
sure gradient curve in such a way, then; since the equation 
armonic may be treated separately and the results added. 


In all previous applications of Fourier series, t 
analytically so that the components could be obtained by integration either analytic- 
ally or numerically; whichever way we used, each component could be obtained to 
any required degree of accuracy. The present situation is different; the pressure 
gradient is given numerically by an experiment so that it is necessary to use a numerical 
technique for determining the components. However, there are experimental errors 
so that only the first few terms of the series can be found accurately; the remaining 
terms, being small, have the same size as the experimental errors. In order to compute 
the coefficients it is convenient to write Equation (1) in the form 


he function being analysed was given 


др ы А 

a. = 4o + У (a, cos nct + b, sin пал), 

ax п=1 
although (1) is the best form to use for solving the equation. The coefficients of the 
first four harmonics for the pressure gradient shown in the preceding figure are given 
below. (The nth harmonic is a, cos nwt + b, sin nox.) The effect of summing these 
four harmonics is shown in the next figure. 


Fourier analysis of the pressure gradient (units mm Hg/cm) 


i а b; 

1 +1.024 +0.240 
2 —0.126 + 1.346 
3 —0.819 +0477 
4 — 0.305 +0.002 


49 = mean value = +0.335 mm Hg/cm 


Sum of the first two harmonies (1) and third harmonic (2) 


The first and second harmonics (1) & (2) 


w 


360 


T—I—T 
200 —7 240 


280 


н ras 
120 — 1605 300 _ 240 


Sum of first three harmonics (1) and fourth harmonic (2) Sum of first four harmonics 
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16.4.2 Oscillatory Flow Along a Cylinder 


In Section 16.2.2 (SAQ 5) we derived the equation of motion for time-dependent flow 
along a cylinder, 


vx te(—x. 0), re(0 К). Q 


where v = p/p and R is the radius of the cylinder. We now consider the special case 
where the pressure gradient is given by 


= LI = Ag", 3) 


where A is a constant. When this solution has: been obtained we can use the results 
of the previous section to determine the flow for a more realistic pressure gradient. 


We expect the solution to be periodic in response to the forcing term (3) so that we may 
write 
wir) = We: 


when this is inserted into Equation (2) we obtain 


ttr tL O<r<R, (4) 
with the boundary conditions 

W(R) =0 
and 

W(r) is bounded as r ~ 0. 
SAQ 8 


By making a suitable transformation, show that Equation (4) may be written in the 
form 


аш ldu 
props hus 
and hence show that the solution to (4) is 
A Jo(Bar/R 

voca s 
where a = R(ojv)*. 
HINT: See W: page 179. 
(Solution on p. 33.) 
SAQ 9 


Using the result of SAQ 8, the full solution to (2) with the pressure gradient (3) and 
the boundary condition w(R, t) = 0 is 


Jo(ar/R)] , 
(r,t) |1 ші pivot, 
(0) mal Јо) 7 (6) 
If z is small Jo(z) can be approximated by the first few terms of its series expansion 
(Equation (40.4) in W: page 179), 


=й = ® E 
Л) -1-5 + 069) 


Using this result show that, in the appropriate limit, the solution (6) agrees with the 
steady state solution obtained in the solution to SAQ 6. 


(Solution on p. 34.) 
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One of the quantities of interest is the total flux or rate of flow, O(a, 1), which is given by 


R pin 
O(a, t) = [| | wr. Or dr 40 
0 0 


1 
= элде = | yW(Ry)dy, 
9 


where we have substituted y — r/R. Note that the flux is given by a complex expression 
Since we are working with a complex pressure gradient. The physical interpretation 
is that the real and imaginary parts of 0 give the fluxes corresponding to pressure 
gradients which are the real and imaginary parts of др/дх, respectively. 


SAQ 10 


(a) Show that the rate of steady flow of a viscous fluid along a cylindrical pipe of 
radius R, and with a constant pressure gradient — A, is 


The suffix s is used to denote the fact that the flow is steady. 


(b) Using the relation 
1 
1 
[ео = 50. 
o č 


which can be deduced from Equation (40.6) in W: page 180. show that 
mR?A [ 2J,(3a) Je 


Qla, t) = 


iwp 7 вәлә) 
_80, ARGON ы 
= el magis |^ (7) 


(Solution on p. 34.) 


In general the terms Ј,( їо) cannot be computed from a simple formula and ап 
involved numerical procedure must be used; in the appendix we give graphs for the 
real and imaginary parts of these functions in the cases n = 0, 1. However, when 
a « lora » 1 there are closed forms which are good approximations to the function. 
For æ « 1 we may use the first few terms of the series expansion 


ЕТТІ БЕЛ 
ER S d eo 4) |. 

H0 = сэ [ % аа 267 
as іп SAQ 9. 


54011 
Show that for small « (low frequencies) 
O(a, t) = Qe", 


ie. when the period is long, the fluid has time to adjust to the slow changes in the 
pressure gradient. 


(Solution on p. 35.) 


If on the other hand х » 1 we may use the asymptotic form of the Bessel function 
(see, for example. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions 
(Dover, 1965), p. 364). 


abbr 


to give, after some manipulation. 


1,5 |а 4п — 1 
Ја) ~ ous ез" expi + 7$ т[ 
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SAQ 12 


Show that for « » 1 


Ola, t) ~ 50. |: - see jq 


io? 


8Q. ui un 
LOTES eilor- n1) 
a 


ie. if the pressure gradient changes rapidly, the fluid motion is always out of phase 
and behind the pressure fluctuations. 


(Solution on p. 35.) 


Having determined Q in the form of Equation (7) it is necessary to evaluate the 
expression. We need to know the typical values of g, some of which are shown in the 
table below. 


Some values of x in different mammalian species 


Species Pulse rate (тіп” 1) Radius (cm) а = К(ш/\у)ї 
Site — root of aorta 
mouse 600-730 0.03 -0.04 1.19- 1.74 
rat 360-520 0.045-0.095 1.38- 3.5 
rabbit 205-220 0.17 3.92- 4.07 
cat 180 0.2 44 
dog 72-125 0.55 -0.6 8.27-10.68 
man 55-72 1.08 -ІЛІ 13.5 -16.7 
ox 43 20 214 
elephant 40-50 447 48-51 
Site — femoral artery 
rabbit 210-360 0.06 14 - 18 
dog 72-180 0.12 -0.15 1.65- 3.25 
man 60-72 02 -035 25-35 


The range of average body weights of these species is from 35 р (mouse) to 
2.0-2.5 х 10% (elephant). The higher range of pulse rate allowed in the 
femoral artery of the rabbit and dog are based on observations in anaesthetized 
animals. 


По = 1, which it is in a large number of cases, there is по simple approximation to the 
Bessel function, and it is necessary to resort to tables. 


Using Equation (7) and the Fourier decomposition of the pressure gradient we can 
find the rate of flow for any given periodic pressure variation. In the first instance 
we notice that the mean value of О(о, t) is zero except when « = 0 (in which case 
it is Q,). so that none of the periodic components of the pressure gradient contributes 
to the overall flow. However, they do contribute to the energy dissipated by the system, 
as we shall now see. 


PDE 164.3 
16.4.3 Dissipation of Energy 


We now consider the energy required to maintain the flow. The energy per unit volume 
is simply the kinetic energy 1pw?, so that the energy in a tube of length | is 


x f pR 
E(t) = i | (#р%»?)2лг ar} dx. 


The energy itself is of no interest; the interesting quantity is the mean rate of change 
of the energy. which measures the average rate at which the heart is doing work on 
this part of the system. Thus we first calculate the rate of change of energy 


dE ШАЙТ; 
ос 9 oF ar Е 
dt тр І Í, w д r dr dx. 


SAQ 13 


Using the equation of motion (Equation (2) in Section 16.4.2) and the equation above 
show that . 


dE dE? ЕФ 
асса ta 
where 
dE" 
d 


= -Q(x. t (px + 1.) — p(x, 0) 


and 


dE? R [ow|? 
a -2ліш [ (=) rdr. 


(Solution on p. 35.) 


The first of these terms represents the rate at which work is being done to maintain 
the pressure difference; if the pressure is periodic in time with zero mean value the 
mean of this quantity is zero. The second term is the rate at which work is being done 
against the viscous forces; it is always negative (unless w = 0). 


In deriving the result in SAQ 13 it is implicitly assumed that the velocity w is real; 
in our calculations w is complex because we chose a complex pressure gradient in 
Equation (3). To make the velocity real we consider the pressure gradient 


= = — Ae" — Ae7i"' = 2[Im(A) sin wt — Re(A) cos wt] 


which is real. Since the equation of motion is linear, the velocity is, for this pressure 
gradient, 


wr, t) = W(r)e'*' + Wire, 


which is also real, so that 


aw\? dw)? ам)? sio dw? 2t (8 
Е -2% dos em 2) е ) 


We consider now the mean value of dE'?)/dt over one cycle, 


ағ” R [ew 
( di )- 2и | (6) ) dr, 


where < ) denotes the time average defined by 


1 T 
qo = 3]. ШОЛ 


T being the period. Since (е = 0, we find, using Equation (8), that 


2 
r dr. 


2 


= к\р 
( dalo E dr 


If we substitute for W using the result of SAQ 8 we find that 
2| z 
$ 


с А 
= —4nl 2. 
( di Шар |748 


ор 
where 2 = іо. On using the result 


1 1 
| ул (uy) y)dy = (ил), (0) — 9707,0, 
9 


и? = v? 


2 pi Е 
NS 
0 


which can be obtained, using Equations (40.5) and (40.6) in W: pages 179 and 180, 
after several lines of manipulation, we find that the mean energy dissipated per unit 


time per unit length is 
2 " 
zT 
R = |. 9 
"n n 


14ЕФ Rt |A£J,(&) 
жа Д 
l dt pat 
(a) By using the recurrence relation (which may be obtained as in SAQ 2 of Unit 14) 


ЕК) 


540 14 


Ja uat = T qa 


with z = č = ia, where « is real, show that Equation (9) reduces to 


Qu me. ARGIN 
Га/77а “Го 


(b) Deduce that in the limiting cases « « 1 and a >» 1 we obtain 


( EN т|А|?К* 
-——)х— а «1, 


ШЕГІ ар 
a Da ЕЕЕ МЕГЕН 
l d/` pet 


(Solution on p. 36.) 
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16.5 COMPARISON WITH EXPERIMENTS 


16.5.1 Sinusoidal Flow 


For our first comparison we consider the simple case where the pressure gradient is 
purely sinusoidal, 


др А ith 
== = asin wt wit 
эх атеа], 


ѕо that the velocity and flow rates аге given by the imaginary parts of Equations (6) 
and (7) in Section 16.4, Iespectively. There are two reasons for making this comparison. 
The first reason is that, until the experiment quoted here had been done, there were 
few experiments measuring the properties of time-dependent viscous flow, and the 
second is that this experiment provided a good check on the apparatus which was 
subsequently to be used for more complicated pressure gradients. 


To compare the velocity result with experiment we need to evaluate the Bessel func- 
tions contained therein; we shall not go into this detail here. The numerical values 


chosen for the experiment ought to be fairly close to the scaled values determined in 
Section 16.3; we choose the following. 


a = 207.7 N m^? = 1.56 mm Hg/m 


T = 6.4 sec 

w = 0.98 rad/sec 
и = 260cp 

R = 2.54 cm 

« = 5.38 


The experimental results and the numerical evaluation of the velocity from the 
theoretical solution are shown in the next figure. In general we see that the agreement 
is good; near the beginning of the cycle the experimental results are ahead of the 
theoretical ones, whilst the converse is the case towards the end of the cycle. This is 
probably due to experimental error. 


Let us now consider the rate of dissipation of energy due to viscous forces. For the 
experiment considered above а = 5.38, and for such large а each Bessel function is 
given by the first term in its asymptotic expansion to within about 5%. Thus, 
with reasonable accuracy, we may use the formula obtained in SAQ 14. Substituting 
the above numbers, we find that 


(2) 2p4 
ae Ya AR where A = 0.5a 
l dt ТАРЫ 


107? watts m^! 


R 


= 4 x 1076 horse power ft^ !, 


which is very small. For the arterial system, using the data given in Section 16.3, 
and assuming a pressure gradient variation of 3mm Нест”! we get a figure of 
about 107* h.p. ft^! ; this is again negligible. These results suggest that there are 
other, better mechanisms for absorbing the energy; one would be the flow through 
small capillaries. 
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Comparison of theoretical and experimental velocity profiles for sinusoidal flow in a rigid pipe. 
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16.5.2 Actual Observations 


In order to compare our theory with experiment, we first measure the pressure gradient 
and perform a Fourier analysis on the result to obtain 


др ы Е 

~ (2 шысы > (a, cos not + b, sin пал) (1) 
х=хо п= 

where N is some number (about 4) determined by the accuracy of the experiment. 

In Section 16.4 we calculated the flow Qla, t) for the complex pressure gradient 

— Ae. It follows by linearity that the rate of flow, with the pressure gradient given 

by Equation (1), is 


00) = 0, + У (a, Ке[О(о /n, 0] + b, Im[Q(v./n, 1). Q) 


п=1 


In the next figure we show а comparison of (2) with actual observations on a dog; 
in this case « = 2.7 and our results are taken from Womersley (see Bibliography). 
Considering the simplicity of our model the agreement is remarkably good. 


Comparison of calculated and observed flow in the femoral artery of the dog. 


Broken curve - observed 
Continuous curve - calculated from observed 
pressure-gradient, x = 2.7 
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16.6 SUMMARY 


Ап attempt has been made to indicate how some of the mathematics developed in 
this course is used in practice. The work described here is necessarily of limited validity; 
the extension that takes account of the elasticity of the wall is made in the paper by 
Womersley. However, as we saw in Section 16.3, there are many other inaccuracies 
inherent in this model. 


It has been shown that mathematics on its own is of limited use in obtaining insight 
into a physical situation and that only by using a variety of techniques do we begin 
io understand a physical system; this simple fact is neglected in some applied 
mathematics texts. 


The methods of attack in the two case studies, Unit 7, Overhead Wires and this unit, 
are typical of the approach towards solving physical problems by use of the applied 
mathematics developed in this course. 
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16.7 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 


In this case we consider the y-component of the equation of motion which is 
T 
p Bp y + viscous force), 


since (v -grad)v = 0 when v(r, t) = (0, 202, 2), 0). For this flow we have shown that 
the viscous forces are given by Equation (6) in Section 16.1, so that the full equation 
of motion becomes 


Solution to SAQ 2 


(a) The solution follows as in SAQ 1 by noting that the pressure is constant, that the 
external forces are zero, and that (г, t) = (0, v(x, t), 0. Thus we have 
9v иде Pv 
0t pox? “ра 
On the boundary x — 0, the fluid velocity is the same as that of the boundary, 
so that 
v(0, t) = u(t). 


(b) It is clear that, because of the oscillatory motion of the boundary, the fluid 
motion at any point will be oscillatory, and because the equation is linear the 
frequency of the oscillation will be that of the boundary. The problem is simplified 
by imposing the complex boundary condition a(t) = Ae’ and taking the real 
part of the solution obtained with this boundary condition. 


Since we expect an oscillatory solution we consider a solution of the form 
d(x, t) = У(х)е'®', 


which gives, on substitution into the equation of motion, 


2 
REC ішУ = 0. 


The general solution to this is 


i 
V(x) = «ep - (5а + әх | + ОЕ (1+ | 


and, since i(x, г) is bounded as х ~ oo, we must have Й = 0; then the boundary 
condition at x = 0 gives a = A, and 


B(x, t) = A exp(—kx) exp Қол — kx) 


+ 
[^] 
Ms Іш 


The real part of this is 


where 


v(x, t) = Ает“ cos(wt — kx), 


which represents a wave travelling into the fluid with speed w/k = (2vcj)* and 
with an exponentially decaying amplitude. 


The force per unit area is given by Equation (3) of Section 16.1 as 
Qv . 
^ Ox} x=0 
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But 


% = kAe ^" (sin(ct — kx) - cos(wt — kx)} 
CX 


= —J/2Ake^* cos(wt — kx + $n), 


so that 


0х 


This force is ош of phase with the motion, i.e. the maximum amplitude of the 
force does not coincide with that of the motion. 


m (5. = —\/2Akp cos(ot + 1л). 
x=0 


Solution to SAQ 3 


(a) (i) The motion is clearly independent of the time so that 00/0; = 0 and Equation 
(2) in Section 16.2.1 reduces to 


v v 

P 

w+- 3-0. 
r r 


The general solution to this is 
B 
v(r) 2 Ar + x 


and the boundary conditions v(a) = 0 and v(b) = bQ give 


Aa? + B = 0, 

Ab? + B = Qb?; 
so that 

Qb? —Qa?b? 

А= р.ғ ра" 
апа 

al) = Qb? в- a? | 

b-a Ë 


(ii) In this case the solution is a little more awkward. The problem is similar 
to SAQ 2(b) and we use the same method to solve it. We impose a boundary 
condition at r = b which is complex, 5(b, г) = АЬе'®', and extract the real 
part of the ensuing solution. As before we suppose that 

B(r, t) = V (rje, 


so that substitution into Equation (2) gives 


шу ye ELE 
н Р fe i 
or 
dV dV 
22 2 2 
2 J7 +15 tE )v-0 
where z = dr, 2* = —iwp/p. This is Bessel's equation of order 1, which 


yields the general solution, 
V(r) = aJ (år) + ВУ, (Ағ). 

The boundary conditions at r = a, b then give 
aJ (4a) + BY,(Aa) = 0 

and 


aJ Ob) + BY,(Ab) = Ab, 
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(b) 


(с) 


so that 


а М0), Ab Q0) 
y Ne 


where 


y = J4(Ab) Y, Qa) — J,(2a)Y, (АБ). 


The solution is thus 
x Ab s 
(r, t) = * {J,(Ar) Y Qa) — У,(2ғ)7 ,(2a)} e. 


The solution for our problem is then the real part of this, which may be 
obtained numerically. 
The moment on the inner cylinder is (per unit length) 
эла? 9v _ 4mpQb?a?- 
or},-, 2 – а? 
The kinetic energy of an elementary ring of fluid of width Ar is ` 
3 mass x (velocity? = 4(2л>Агр) + v? 
= пру? Nr. 


The total kinetic energy is thus 
b 
E= np | rv? dr. 


The rate of change is 


dE ар 
— = е rv? dr 


dt dt 
= wf 07 
= 2лр ШЕ 
Substituting for dv/dt from the equation of motion we obtain 
= = anu [ |» a + v - La di 


We integrate (by parts) the first term in the integral and obtain 


dE , I-II + ac | Qv о? di 
а 2" Hal, 22” a} ы "Or т] 


=2 Pd > [ |ы 
= nu "|, ли i АЕ „| 


We consider the first term. The force per unit area of the boundary is и Qv/Cr, so 
that the rate of working of the boundary is д(ь Ov/Or),-, per unit area, where 
0 = a or b. Thus the rate of working on a cylinder of length 115 2л/шге Cv/r. 
This term thus represents the energy given to the fluid by a moving boundary. 
If the boundaries are stationary, v(a, t) = v(b, t) = 0 and no energy is lost or 
gained. 


The second term is simply the energy lost owing to viscous forces in the liquid; 
the integrand is always positive and so this term is always negative. Thus the 
energy is always decreasing, and the fluid will ultimately come to rest; the energy 
is not annihilated but is converted into random motion, which means that the 
fluid gets heated. 
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Solution to SAQ 4 
The viscous force is obtained from the formula 


force = и x area x velocity gradient 


at a radius 7. On the inner wall F = r — Ar and on the outer wall ? = r + Ar, whence 
the result quoted. The only other force acting on this element of fluid is that due to 
the pressure at either end, i.e. 


pressure x area, 
where the area is given by 2nr x 2Ar. Thus the total force in the positive x-direction is 


ôw r+ Ar 
— 4лгА [р(х + 1) — p(x)] + 2nul "ж . 


= Ағ 


Solution to SAQ 5 
Newton's Second Law of Motion is 


mass x acceleration — force, 


ie. 
" Ir*ár 
алара 2* = -4лғАғ(р(х + 1) — p(x)) + 2nul бш * 
д or РА 
First we divide by 4zrlpAr and let | ~ 0; then, since 
lig PG 0 = PON др 
1-0 1 дх 
we obtain 
aw — lóp H ,v prar 
atop ax 2ғрАғ| Or |, 
Now proceeding to the limit as Ar ~ 0 and noting that 
| 5] | ғы) 
А ‘Or or, ô | dw 
1 r*ár rar] А ОАТ 
am, 2Ar or | z 
we obtain 
д» _ 1др p д» 12}. 
ді рдх plor? тд 
Solution to SAQ 6 
For a steady flow w is independent of the time; it is given that др/дх = — Ap, so that 


we need to solve 


aw 1 Ow Ap 


а ror Т зау 
or 
d{ dw 
т | x) = ~kr 
The solution is simply, 
г? 
w(r) = — — + clnr 4 d, 


4 
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where c and d are arbitrary constants. Now w must be finite at the origin, so that 
c = 0; also w(R) = 0 so that d = 4kR?. Hence 


k 2 2 Ap 
(r) = — =p) = 2P R2 3 
w(r) 2% ыты ғ). 
Note that the equation for steady flow is just the equation 
Vw = —k 


quoted in Section 3.2.2 of Unit 3, Elliptic and Parabolic Equations, with circular 
symmetry. 


Solution to SAQ 7 
The equation of motion is 
де 1др е E 1 x] 


Ot рдх ôr? ЫЗ T Or 
where v = р/р. On putting r = о, t = fi, it is easy to see that, if (7,7) = war; ft), 
aw 1{, др „ [22% 10$ 
= +7 ---- , 
| ' Е F a 


with 


is the equation describing the fluid flow of the model. The latter equation describes 
flow with kinematic viscosity given by 


j=— 
a? 


and pressure gradient | др/дх. Thus if the pressure gradient is oscillatory with 
period T in the artery it must have period T/f in the model, and its amplitude must 
be a factor fj times the arterial amplitude to yield the same velocity profile. 


Solution to SAQ 8 
Consider the substitution z = (— ic/v)r as in Unit 14 so that Equation (4) becomes 


аи 1 du A 


22 zd "T iop 


where W(r) = u( / —iw/v r). The general solution to this, which is bounded as z ~ 0, is 


A 
u(z) — ішр + cJo(z), 


so that 


Wr) = 5 + cJ (оу) — since(- D$ = (P) = й 


ES E + cJo(itar/R) 
iwp 


where in this last equation we have set 
a = R(w/v)}. 
The boundary condition, W(R) = 0, then gives 
A UEM 
E р 29 — Е 
Win Zl Jo (ido) 
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Solution to SAQ 9 


Using the expansion we have 


A iar? io? 
«ep bre e] 


where the last expression is obtained by using the binomial expansion and ignoring 
terms of order a^. Expanding the brackets and again ignoring all terms of order a* 
we obtain 


Aa r? 
wisis pg) 


so that 
wr, t) с fr- Pje, 


Now one way in which о сап be small is if c « 1; thus, if the period of the oscillations 
is large, the fluid motion is approximately in phase with the pressure changes, i.e. 
the fluid has time to adapt to the changes. This would be expected on common-sense 
grounds. 


In the limit as о ~ 0 we have the steady state problem, and our solution agrees with 
that obtained in SAQ 6, with А = Ap. 


Solution to SAQ 10 
(a) For steady flow we use the result of SAQ 6, with А = Ap, 


w(r) = fr = n) 
The rate of flow is now 


0, = 2n Г w(r)r dr 
0 


= (R?r — r°) dr 
ЛАҚ 
8и 
This is of course the résult which you obtained іп SAQ 14 of Unit 3. 
(b) The total flow is 


00, _ 21А giat f b E Д dy 


dep 0 Jolita) 
_ TRA ja 2J (8а) ET 
ЕСТЕ Bad (ita) 
on using the given relation. Since, from (a), 


nAR* 2 пЕ?А іо? 
80 шр 8 


0,- 


the result follows. 
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Solution to SAQ 11 


Consider the expression 


which we obtain on using the small а expansion. Using the Binomial Theorem this 
is equal to 


ім 


‚Юю? ig? io? 
1 ! + ehi- =| + 0 = z + 008“), 


and the result follows using the formula obtained in SAQ 10(b). 


Alternatively, the result may be obtained directly by using the result of SAQ 9. 


Solution to SAQ 12 


Using the asymptotic form quoted we find that for a » 1 


so that, from Equation (7) in Section 16.4.2 we have 


80 зе") 50 ы 24 
06) v2 тұла |е“ = Bees м), 


80. , ЕТ 
Qla, t) ~ 80. eead since j = e™2, 
a 


Solution to SAQ 13 


Substituting the equation of motion into the expression for dE/dt we can separate 
out two terms, 


(1) xl 
“= -2 f Г ДЕ Pares 
d 


хк д 
-2x |" w(r)r dr [ 554 


— Q(a. (р(х + 1,0) — p(x.) 


апа 


dt or 
Е м ya ду “ar 
= zz w(r) Or 
since w(r) is independent of x. Integrating the first term of this last integral by parts 
we have 
dE? ow |^ [ wo “© j- EIL ] 
-2 "Ww 
di 2nyl 5| rw rs or 
R [Ow 
—Ó --|ра 
шығы ji { 5; " 


since w(R, t) = 0. 


(2) xl 
ae = anu f [ o - +129 
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Solution to SAQ 14 
(a) Putting n = 1 in the recurrence relation, we have 


2 


(ш) = =J (2) — Jolz) 


so that 


922 _ Jol) 
UD C UO 


and since č = ita, č? = — іа? is purely imaginary so that 
р gl 


AG) KS 1 
= —Re М 
ке 9 (5 - 73) d 5. sy 


1 
Re [Б] =® ee) =ъ jor Re(®). 
Thus 


RAG 2) ШЗ 
EXC |249. dud Етті 


and the result follows. 


(b) Ifa « 1 we have 
EQ cea - ae 
Jol) 1-40 


—Aio(1 + dic?) (1 — 1092) 


using the first two terms of the series expansion 


R 


2 
ы 
ҮЧ 


= —tie? — та“. 
Thus 
(um. ADR. at n|APR*. 
at 16) | 4u 


If a > 1 we have, as in the solution to SAQ 12, 


GI RN T NE ey Di) o а 
ЕТЕ = Re(ce'"^) = Re(— ita) = —«Re =r ia - я 
Thus 
14Е® m An|AP R*. [2 4n|APR* 
гау pot — 23 цї?2} 
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16.8 APPENDIX 


Bessel Functions 


Here we try to give you an idea of how the Bessel functions J9(i#x) and J (х) behave. 
These functions have complex values and we examine their real and imaginary parts 
separately. 


In the first two figures, we have plotted the real and imaginary parts of the three 
functions 
Дх) = Jo(ix)e7 92, 


„2 
"E 
4 
which approximates / by taking the first two terms in the expansion of Jo, and 


h(x) = S expi a j 


57% 


obtained using the asymptotic form preceding SAQ 12. 


0.7 


0.64 Im Joli? ?x) e7 7^? Re Joli? 2 )e7 7? 


0.54 


We see that the power series expansion is good up to about x — 2, after which it 
becomes progressively worse. For x > 2 the asymptotic form for the imaginary 
part is indistinguishable from the correct values, and for the real part the same is 
true for x > 4. This is shown quite clearly in the following table. 


x Re fo(x) Re hox) Im f(x) Im ho(x) 
10 0.485 0.379 0.123 0.123 
20 0.183 0.147 0.236 0.241 
3.0 — 0.0265 — 0.0362 0.232 0.227 
40 —0.152 —0.152 0.136 0.129 
5.0 — 0.182 —0.178 0.0034 — 0.0002 
6.0 —0.127 —0.124 —0.105 — 0.106 


Notice that there is a region around x = 3 where both 20 and л, are poor approxi- 
mations. It is typical of these types of approximations that there is a region where 


neither works very well. 
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In the next pair of figures we show 
fio) = J lixe? 
together with the form derived from the series expansion 
x 


(x)= ire - ons 
gx WA gti zle 


and the asymptotic form 
I [x Зл 
hy(x) = aux expi JA S 3 


0.24 ReJ,(i? tne? 0.2 


Im J, (33 2 x)e7 ^? 


The comments for these figures are much the same as for the previous pair, but notice 
that the asymptotic expansion is not very accurate until x == 5 for the real part and 
x = 7 for the imaginary part: this is shown clearly in the table below. 


х Re f,(x) Re h,(x) Im / (х) Im л, (х) 
10 —0.195 —0.123 0.152 0.379 
20 —0.242 —0.241 0.0739 0.147 
30 —0.208 —0227 — 0.0584 —0.0362 
40 -O111 —0.129 —0.152 —0.152 
5.0 0.0105 0.0002 —0.169 — 0.178 
6.0 0.107 0.106 —0.113 —0.124 
70 0.144 0.149 —0.0164 — 0.0233 
80 0.114 0.120 0.0757 0.0739 


9.0 0.0357 0.0408 0.124 0.127 
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